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ABSTRACT 


In  this  study,  a  numerical  weather  prediction  model  is 
developed  and  tested  on  some  selected  synoptic  situations 
over  Western  Canada.  A  four-level  model  of  the  atmosphere 
is  formulated  in  which  a  prediction  is  obtained  by  numerical 
solution  of  the  quasi- geostrophic  omega  and  vorticity 
equations.  The  horizontal  grid  (with  a  grid  point  spacing 
of  200  km)  encompasses  Western  Canada  and  parts  of  the 
Northwest  Territories,  the  Yukon,  Alaska,  the  Eastern 
Pacific  Ocean  and  the  northwestern  United  States.  The 
Adams-Bashf orth  time  stepping  scheme  with  a  time  increment 
of  1/2  hour  is  employed  to  achieve  the  forecast.  Three 
synoptic  situations  were  chosen  in  which  the  western 
Cordillera  of  North  America  appeared  to  have  some  influence 
on  the  development  and  movement  of  cyclones  in  Western 
Canada.  In  data  sparse  areas  of  the  grid,  bogus  data  are 
generated  in  an  attempt  to  provide  an  accurate  analysis  for 
initialization  of  the  forecasts. 

The  study  focuses  on  the  importance  of  different 
formulations  of  the  lower  boundary  condition  in  the 
numerical  model.  Terrain  elevations  at  each  grid  point  are 
obtained  from  topographic  charts,  while  frictional  drag 
coefficients  are  interpolated  to  the  grid  from  Cressman’s 
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(1960)  data.  Several  lower  boundary  conditions  which  have 
been  employed  in  other  numerical  models  are  compared  with 
each  other.  The  relative  accuracies  of  the  12  hour 
forecasts  are  obtained  through  the  calculation  of  root-mean- 
square  errors  between  the  forecast  and  the  observed 

i 

geopotential  height  fields. 

The  testing  shows  that  f rictionally-induced  vertical 
motion  at  the  lower  boundary  is  as  important  as  orographic 
lift  and  subsidence  in  accurate  prediction  of  lee 
cyclogenesis.  The  best  forecasts  were  achieved  if  a 
slightly  smoothed  terrain  was  employed  and  the  vertical 
velocity  at  the  lower  boundary  was  assumed  to  originate  at 
the  terrain  surface  rather  than  at  the  1000  mb  surface.  If 
the  850  mb  geostrophic  wind  rather  than  a  geostrophic  wind 
at  terrain  height  was  employed  in  calculations  of  frictional 
and  orographic  components  of  the  vertical  velocity  in  the 
planetary  boundary  layer,  a  more  accurate  forecast  was 
achieved . 

In  conclusion,  the  four-level  guasi-geostrophic  model 
with  a  suitable  lower  boundary  condition  shows  skill  in 
short  term  weather  prediction  over  Western  Canada. 


vi 


' 


ACKNOWLEDGEMENTS 


I  wish  to  express  ray  thanks  to  the  many  people  and 
organizations  who  made  completion  of  this  study  possible. 

For  his  helpful  guidance  in  the  latter  terra  of  this 
study,  I  wish  to  thank  Dr.  Edward  P.  Lozowski.  A  special 
note  of  thanks  is  also  due  Dr.  Madhav  L.  Khandekar  who 
provided  the  initial  impetus  in  the  formulation  of  the 
problem.  I  also  thank  Dr.  Keith  D.  Hage  and  Dr.  Stan  Cabay 
who  served  on  ray  examining  committee  in  addition  to 
Dr.  Lozowski  who  acted  as  chairman. 

I  gratefully  acknowledge  the  Division  of  Meteorology 
for  providing  the  funding  for  the  extensive  computing  which 
this  study  entailed.  The  helpful  advice  offered  by  the 
staff  and  many  of  my  fellow  students  in  the  Division  of 
Meteorology  is  also  sincerely  appreciated.  The  excellent 
photo-reduc tion  of  the  many  diagrams  was  the  responsibility 
of  Mr.  J.  Chesterraan. 

This  study  was  conducted  while  I  was  on  Educational 
Leave  from  the  Atmospheric  Environment  Service,  Environment 
Canada. 


Vll 


■ 


) 


TABLE  OF  CONTENTS 


Page 

DEDICATION  iv 

ABSTRACT  v 

ACKNOWLEDGEMENT  vii 

TABLE  OF  CONTENTS  viii 

LIST  OF  TABLES  xi 

LIST  OF  FIGURES  xii 

CHAPTER 

I  INTRODUCTION  1 

1.1  Preliminary  Comments  1 

1.2  Studies  of  Lee  Cyclogenesis  and 

Numerical  Modeling  2 

1.3  Outline  of  This  Study  5 

II  THE  MODEL  6 

2.1  The  Prediction  Equation  6 

2.2  The  Diagnostic  Omega  Equation  8 

2.3  The  Horizontal  Grid  9 

2.4  The  Vertical  Representation  11 

2.5  The  Lower  Boundary  Condition  13 

III  DATA  ACQUISITION  AND  ANALYSIS  20 

3.1  Upper  Air  Data  20 

3.2  Surface  Data  22 


Vlll 


Page 


3.3 

The  Objective  Analysis 

26 

3.4 

The  Static  Stability  Parameter 

28 

NUMERICAL  PROCEDURES 

31 

4.1 

Horizontal  Finite-Difference  Operators 

31 

4.2 

The  Vertical  Derivatives 

36 

4.3 

The  Time  Step 

38 

4.4 

Solution  of  the  Differential  Equations 

43 

4.5 

The  Forecast 

50 

RESULTS  AND  DISCUSSION 

52 

5.1 

Preliminary  Comments 

52 

5.2 

Testing  of  Numerical  Procedures 

53 

5.3 

Case  1:  January  5,  1972 

56 

5.4 

Comparison  of  Several  Forecast  Runs: 

Case  1 

63 

5.4.1 

Smoothing  of  Topography:  Runs  1A,  IB, 

1C 

63 

5.4.2 

Friction  and  Topography:  Runs  2A,  2B, 

2C 

73 

5.4.3 

Small  Frictional  Vertical  Velocity: 

Run  3 

84 

5.4.4 

The  850  mb  Wind  in  PBL  Computations: 

Run  4 

87 

5.5 

Case  2:  March  5,  1972 

90 

5.6 

Case  3:  April  27,  1972 

97 

5.7 

Systematic  Errors  in  the  Model 

104 

ix 


VI 


SUMMARY  AND  CONCLUSIONS 


109 


BIBLIOGRAPHY 

113 

APPENDIX  A 

The  Aitken  Polynomial 

117 

APPENDIX  B 

Spatial  Smoothing 

121 

x 


LIST  OF  TABLES 


Table 

1 

Surface  regimes  and 

drag  coefficients 

Page 

25 

2 

The  static 

stability 

parameter 

30 

3 

Statistics 

of  Case  1 

,  Run  1 

67 

4 

Statistics 

of  Case  1 

,  Run  2 

74 

5 

Statistics 

of  Case  1 

,  Run  3 

86 

6 

Statistics 

of  Case  1 

,  Run  4 

89 

7 

Statistics 

of  Case  2 

96 

8 

Statistics 

of  Case  3 

102 

xi 


LIST  OF  FIGURES 


Figure  Page 

1  The  horizontal  grid  area  10 

2  The  vertical  structure  of  the  model 

atmosphere  12 

3  Relationship  between  the  geostrophic 

wind  and  the  wind  in  the  PBL  17 

4  Pressure  at  terrain  height  23 

5  Cressman  drag  coefficients  24 

6  Horizontal  grid  point  structure  32 

7  Comparison  of  centered-dif f erence  and 

Adams-Bashf orth  time  step  procedures  44 

8  Initial  analysis  at  850  and  500  mb  for 

Case  1  57 

9  Initial  surface  analysis  for  Case  1  58 

10  Verification  analysis  at  850  and  500 

mb  for  Case  1  60 

11  Verification  surface  analysis  for 

Case  1  62 

12  Smoothed  pressure  at  terrain  height  65 

13  12  hour  forecast  at  850  mb  for  Funs  1A 

and  IB  66 

14  Same  as  Figure  13  but  at  500  mb  68 

15  Same  as  Figure  13  but  at  the  surface  70 

16  12  hour  forecast  at  850  mb  and  surface 

for  Run  1C  72 

17  Same  as  Figure  16  but  for  Run  2A  75 


Xll 


Figure 


Page 


18 

Initial  vertical  velocity  fields  at 
the  surface  and  at  600  mb  for  Run  2B 

77 

19 

12  hour  forecast  at  850  mb  and  the 
surface  for  Run  2B 

78 

20 

Same  as  Figure  18  but  for  Run  2C 

80 

21 

Same  as  Figure  19  but  for  Run  2C 

82 

22 

12  hour  forecast  at  850  mb  and  at  the 
surface  for  Run  3 

85 

23 

Same  as  Figure  22  but  for  Run  4 

88 

24 

Initial  analysis  at  850  and  500  mb  for 
Case  2 

91 

25 

Verification  analysis  at  850  and  500 
mb  for  Case  2 

93 

26 

850  mb  forecasts  for  Case  2 

94 

27 

500  mb  forecast  for  Case  2,  Run  1 

96 

28 

Initial  analysis  at  850  and  500  mb  for 
Case  3 

98 

29 

Verification  analysis  at  850  and  500 
mb  for  Case  3 

99 

30 

850  mb  forecasts  for  Case  3 

101 

31 

500  mb  forecast  for  Case  3,  Run  1 

102 

1 


CHAPTER  I 


INTRODUCTION 


1 . 1  Preliminary  Comments 

Numerical  prediction  of  meteorological  systems  on  the 
synoptic  scale  has  achieved  a  high  degree  of  sophistication 
over  the  past  two  decades,  and  many  of  the  physical 
processes  which  govern  atmospheric  motions  have  been 
included  in  the  numerical  models.  Only  limited  success  has 
been  achieved  in  the  simulation  of  the  cyclogenesis  observed 
in  the  lee  of  major  mountain  ranges,  although  most  models 
are  able  at  least  to  predict  the  gross  characteristics  of 
this  phenomenon.  Since  advances  in  computer  technology  have 
allowed  the  direct  solution  of  the  primitive  equations  of 
motion,  most  attention  has  recently  focused  on  primitive 
equation  models.  However  such  models,  while  logically  quite 
straight-forward,  tend  to  be  difficult  to  implement  because 
of  their  sensitivity  to  the  initialization  of  the  data.  As 


2 


well  they  require  a  large  amount  of  computing  to  achieve  a 
forecast. 

A  numerical  model  based  on  the  system  of  filtered 
geostrophic  equations  is  somewhat  less  stringent  in  its 
requirements  because  gravity  oscillations  present  in  the 
primitive  equations  models  have  no  solution  in  the 
geostrophic  formulation.  Although  the  geostrophic 
assumption  is  a  powerful  constraint  between  the  horizontal 
wind  and  the  pressure  gradient,  especially  near  the  surface 
of  the  earth,  successful  duplication  of  synoptic  scale 
motion  using  the  geostrophically  filtered  system  of 
equations  has  been  achieved  for  many  years.  In  this  study, 
a  geostrophic  model  with  four  levels  in  the  vertical  is 
developed  and  tested  on  an  area  encompassing  Western  Canada, 
including  parts  of  the  Northwest  Territories,  the  Yukon, 
Alaska,  the  Eastern  Pacific  Ocean  as  well  as  the 
northwestern  Onited  States.  The  purpose  of  this  study  is  to 
assess  the  influence  of  the  lower  boundary  condition  of  the 
numerical  model  in  the  prediction  of  weather  systems  in 
Western  Canada,  especially  in  relation  to  the  simulation  of 
the  lee  cyclogenesis  process. 

1 . 2  Studies  of  Lee  Cyclogenesis  and  Numerical  Modeling 

Since  it  is  through  the  vertical,  velocity  that 
atmospheric  potential  energy  is  converted  into  kinetic 
energy  (Lorenz,  1955),  many  diagnostic  studies  have  examined 
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the  vertical  velocity  or  omega  fields  ( co  =  )  at  the 

d.  t 

surface  of  the  earth  as  well  as  at  levels  well  above 
orographic  features.  In  an  early  study,  Newton  (1956) 
calculated  vertical  velocities  using  the  adiabatic 
assumption.  He  indicated  that  the  horizontal  variation  of 
frictional  drag  at  the  surface  may  be  of  primary  importance 
in  the  initial  low-level  lee  cyclogenesis. 


Diagnostic  studies  using  the  g uasi-geostrophic  omega 
equation  have  been  employed  to  study  the  role  of  the 
vertical  velocity  field  in  the  mechanism  of  circulation 
change  in  the  lower  to  middle  troposphere.  Haltiner,  Clarke 
and  Lawniczak  (1963)  and  Stuart  (1964)  utilized  omega 
equations  to  evaluate  the  co  fields,  and  to  estimate  the 
relative  magnitudes  of  terms  in  the  vorticity  equation 
responsible  for  the  development  of  baroclinic  waves.  In  a 
later  study,  Schram  (1974)  showed  that  with  certain  lower 
boundary  conditions,  areas  of  vorticity  production  were  well 
correlated  with  favored  areas  of  lee  cyclogenesis  in  Western 
Canada. 


Actual  numerical  prediction  using  the  quasi-geostrophic 
filtered  equations  has  been  extensively  developed,  but  few 
studies  have  focused  on  the  lower  boundary  condition  and 
especially  on  the  lee  cyclogenesis  problem.  Many  models 
formulated  for  prediction  on  limited  area  grids  have 
included  crude  pararaete rizations  of  the  processes  within  the 
atmospheric  boundary  layer  at  the  surface  of  the  earth.  The 
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early  two-level  model  of  Thompson  and  Gates  (1956)  included 
only  the  large  scale  ascent  of  air  over  orographic  features. 
The  incorporation  of  a  representative  lower  boundary 
condition  in  numerical  models  was  discussed  by  Sawyer 
(1959).  His  derivation  of  frictional  effects  was  based  on 
the  theory  of  the  flow  in  the  planetary  boundary  layer  as 
well  as  on  synoptic  experience.  The  two-level  model  of 
Greystone  (1962)  was  tested  on  fifty  different  synoptic 
cases  over  the  north  Atlantic  and  Europe.  Only  a  small 
change  in  forecast  error  was  noted  when  forecasts  employing 
a  sophisticated  lower  boundary  were  compared  with  those 
employing  no  orographic  or  frictional  influence. 

Most  filtered  equation  models  have  made  use  of  the  drag 
coefficients  tabulated  by  Cressman  (1960)  for  the  entire 
northern  hemisphere.  As  well,  the  highly  smoothed  orography 
of  the  northern  hemisphere  provided  by  Berkofsky  and  Bertoni 
(1955)  is  generally  employed  in  prediction  models  such  as 
those  of  Cressman  (1963)  and  Danard  (1966).  More  recent 
work  with  the  primitive  equations  has  allowed  the 
incorporation  of  more  sophisticated  boundary  layer  models 
into  the  large  scale  predictions.  The  primitive  equation 
model  employed  by  Egger  (1972,  1974)  was  used  to  simulate 
the  flow  in  the  vicinity  of  the  high  mountains  of  Europe  and 
North  America,  while  Danard  (1971)  studied  the  formulation 
of  frictional  effects  in  the  boundary  layer  of  his  primitive 
equation  model. 
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1 . 3  Outline  of  This  Study 

In  this  study,  a  four-level  model  is  developed 
employing  the  quasi-geostrophic  system  of  equations  to 
obtain  short  range  (12  hour)  prediction  of  weather  systems 
over  Western  Canada.  Testing  of  the  model  was  conducted  on 
three  cases  in  which  some  degree  of  lee  cyclogenesis  was 
observed  to  occur  within  the  short  prediction  period. 

The  emphasis  in  the  study  is  on  the  formulation  of  a 
lower  boundary  condition  for  the  numerical  model  which 
permits  the  most  accurate  prediction  of  the  motion  and 
development  of  weather  systems  in  the  vicinity  of  the 
Canadian  Rocky  Mountains.  Six  different  parameterizations 
of  the  lower  boundary  condition  are  tested.  Quantitative 
comparisons  of  the  forecasts  are  obtained  through 
calculation  of  the  root-mean-square  errors  between  forecast 
and  observed  (or  verification)  geopotential  fields. 

In  Chapter  II  the  quasi-geostrophic  system  of 
prediction  equations  and  the  details  of  the  model  atmosphere 
are  briefly  described.  Chapter  III  describes  the  data 
extraction  and  analysis  procedure,  while  the  numerical 
techniques  employed  in  the  model  are  presented  in  Chapter 
IV.  The  results  of  the  testing  of  the  three  cases  are 
detailed  in  Chapter  V,  and  in  Chapter  VI  a  brief  summary  and 
the  conclusions  of  the  study  are  presented. 
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CHAPTER  II 


THE  MODEL 


2. 1  The  Prediction  Equation 

The  prediction  equation  used  in  this  study  is  an 
approximation  to  the  complete  vorticity  equation  which,  when 
written  in  a  coordinate  system  with  pressure  p  along  the 
vertical  axis  k ,  is: 


+  V*  V(C+f )~ +  x -  O  c2*1) 


In  this  equation,  C”  is  the  vertical  component  of  the 
relative  vorticity,  "f  the  Coriolis  parameter,  \/  the 
horizontal  wind,£->  =  the  vertical  velocity  in  pressure 

coordinates,  and  V  the  horizontal  gradient  operator. 


While  studies  by  Stuart  (1964)  and  Schram  (1974)  have  shown 
that  the  last  term  on  the  left-hand-side  of  equation  (2.1) 
may  on  occasion  be  of  some  importance,  it  is  generally 
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ignored  in  numerical  prediction  models.  The  analysis  of 
Wiin~ Nielsen  (1959)  has  shown  that  if  this  twisting  term  is 
ignored,  the  vertical  advection  term ,  U)  ,  should  also  be 
ignored  since  the  two  tend  to  compensate  each  other.  As 
well,  he  indicates  that  if  the  divergence  term(C  +  -f)  ,  is 

retained,  the  absolute  vorticity,  S  +  f  ,  should  be  replaced 
by  a  constant  in  this  term  only.  All  of  these  modifications 
are  necessary  because  the  integral  of  the  complete  vorticity 
equation  over  a  closed  surface  is  zero.  Hence  any 
simplification  to  the  vorticity  equation  should  attempt  to 
retain  this  property. 


With  these  simplifications,  the  resulting  quasi- 
geostrophic  vorticity  equation  is: 

=  -Va  •  V(C3  +  f)  +  f  fyf  (2. 

is  the  average  value  of  the  Coriolis  parameter 
region  of  interest.  The  horizontal  wind  has  been 
by  the  non-diver  gent  geostrophic  wind: 


where  "f 
over  the 
replaced 


Vj  =  Vf  I  xVZ 


(2.3) 


and  the  relative  vorticity  by  its  geostrophic  counterpart: 


Cj  =  3/f  v2z 


is  the  height  of  an  isobaric  surface 


In  these  equations. 
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in  geopotential  meters  (  gpm  )  and  g  is  the  gravitational 
acceleration  (assumed  constant)  ,  Using  the  standard 
definitions  of  the  Laplacian  and  Jacobian  operators,  the 
quasi-geostrophic  vorticity  equation  may  now  be  written: 

A(vzz)=v"( jp  (2-5> 


2 • 2  The  Diagnostic  Omega  Equation 


In  order  to  derive  a  diagnostic  equation  for  omega,  use 
is  made  of  the  quasi~geostrophic  thermodynamic  equation  in 
the  pressure  coordinate  system: 


ewf  © 

In  this  equation,  =  — 0  is  a  static  stability 

parameter  where  <K  is  the  specific  volume  of  the  air  and  0 
the  potential  temperature.  The  static  stability  is 
permitted  to  be  a  function  of  pressure  only  in  this  version 
of  the  thermodynamic  equation,  for  if  o'  were  spatially 
variable,  it  may  be  shown  (Haltiner  1971,  pp  73*74)  that 
global  integral  constraints  on  the  thermodynamic  equation 
would  be  violated.  With  the  generally  stable  atmospheric 
stratifications  of  winter  months  in  high  latitudes,  o'  is  a 
positive  number,  and  is  evaluated  as  indicated  in  Chapter 


IV, 
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The  time-derivatives  in  equations  (2.5)  and  (2.6)  may 
be  eliminated  by  operating  on  the  former  with  and  on 

the  latter  with  ^7  and  subtracting,  yielding  the  quasi- 
geostrophic  omega  equation: 

-  aj?  J(Z,C3-f)-^J(Z,§|) 

From  this  equation  it  may  be  observed  which  atmospheric 
conditions  tend  to  give  rise  to  negative  CO  (ascent) .  Hith 

is  positive  and 
thus  CO  is  negative.  On  the  other  hand  subsidence  will  be 
associated  with  cold  advection.  If  vorticity  advection 
increases  with  height,  then  T(z,c3+f)  is  positive  and 
ascent  will  be  initiated.  Thus  in  a  situation  with  a  short 
wave  propagating  in  a  westerly  flow,  an  area  of  ascent  will 
be  located  between  the  trough  and  the  downstream  ridge. 

dl 

Equations  (2.5)  and  (2.7)  are  a  closed  set  for  co  and  ^ 
and  are  solved  sequentially  using  numerical  techniques.  The 
approximation  of  these  equations  by  finite-difference 
analogues  and  the  method  of  solution  will  be  described  in 
Chapter  IV. 

2 • 3  The  Horizontal  Grid 

In  the  horizontal,  atmospheric  parameters  are 
represented  on  a  grid  of  points  over  the  area  indicated  in 
Figure  1.  The  grid  is  so  oriented  that  the  prime  region  of 


warm  advection  into  an  area 


J(zj 


Sp 


' 
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Figure  1.  The  areal  extent  of  the  horiz  ntal  grid.  Grid 
point  12,14  lies  to  the  west  of  Edmonton,  Alberta.  The 
scale  indicated  is  the  grid  point  spacing  true  at  60°  N. 
The  area  outlined  is  that  for  which  root-mean-error s  were 
calculated . 
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interest  is  located  near  the  center  of  the  area,  with  the 
Rocky  Mountain  chain  essentially  bisecting  the  area.  The 
grid  itself  is  a  matrix  of  21  rows  and  23  columns  and  the 
grid  element  is  a  square  on  a  polar  stereographic 
projection.  The  grid  interval  is  200  kilometers  (km),  true 
at  latitude  60°  north  while  elsewhere  on  the  grid,  the 
spacing  dj^  is  given  by: 


(2.8) 


where  D  is  200  km  and  m  the  map  scale  factor  defined  by: 


(  +  S  I N  GO 

I  +  SIN  4>(\ 


(2.9) 


In  this  equation 


is  the  latitude  of  grid  point  i,j. 


2 • 4  The  Vertical  Representation 

In  the  vertical,  the  model  atmosphere  has  five  layers 
separated  by  the  four  primary  isobaric  surfaces:  300,  500, 
700  and  850  millibars  (mb)  which  coincide  with  the  standard 
levels  at  which  data  are  measured  in  the  real  atmosphere. 
The  input  of  data  is  at  the  four  primary  levels  and  the 
omega  equation  is  solved  at  the  intermediate  levels,  400, 
600  and  775  mb,  subject  to  the  vertical  boundary  conditions 
that  CO  =0  at  200  mb  and  CO  =  COg  at  the  lower  boundary. 

A  schematic  vertical  cross-section  of  the  model  atmosphere 
is  illustrated  in  Figure  2. 
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Figure  2.  The  vertical  structure  of  the  model  atmosphere 
using  pressure  as  the  vertical  coordinate.  The  vorticity 
equation  is  solved  at  levels  labelled  Z,  and  the  omega 
equation  at  levels  labelled  to  „  The  vertical  velocity  at 
the  lower  boundary  originates  either  at  1000  mb  or  at  the 
terrain  pressure  surface  p  . 
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2 •  5  The  Lower  Boundary  Condition 

The  lower  boundary  in  the  model  is  the  terrain  surface, 
and  the  lower  boundary  condition  is  therefore  the 
specification  of  omega  at  this  surface.  The  vertical 
velocity  has  two  components  in  the  model  which  are 
determined  respectively  by:  a)  the  earth’s  topography,  and 
b)  the  friction  in  the  surface  boundary  layer. 

a)  The  Orographic  Vertical  Velocity 

Haltiner  (1971)  gives  an  expression  for  the  orographic 
component  of  the  vertical  velocity  in  the  form: 

CJt  ^  3  T  Vt*VHtM  (2.10) 

where  H+  is  the  height  of  the  terrain  above  sea  level  and 
ft  and  Vt  are  respectively  the  atmospheric  density  and 
the  horizontal  wind  velocity  at  terrain  height.  Converting 
to  pressure  coordinates,  this  expression  reduces  to: 

cot  =  V-t '  V  /\  (2-11) 

where  p.t  is  the  pressure  at  terrain  height.  In  the  model, 
Pt  is  evaluated  to  be  the  pressure  at  height  H t  in  the 
U.  S.  Standard  Atmosphere.  In  principle,  it  should  be 
possible  to  estimate  terrain  pressure  from  actual 
atmospheric  conditions  aloft,  but  such  a  refinement  was  not 
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believed  to  be  necessary,  since  the  gradient  of  terrain 
pressure  rather  than  its  absolute  value  is  of  importance  in 
equation  (2.11).  In  the  vicinity  of  a  mountain,  the 
synoptic  pressure  gradient  will  be  very  much  smaller  than 
the  variation  of  pressure  with  height.  Thus  the  additional 
computational  complications  required  to  forecast  terrain 
pressures  are  not  justified  by  the  very  small  percentage 
changes  anticipated  in  .  The  details  regarding  the 

determination  of  Ht  and  Vt  will  be  discussed  in  Chapters 
III  and  IV. 

b)  The  Frictional  Vertical  Velocity 

The  f rictionally^induced  vertical  velocity  at  the  top 
of  the  boundary  layer  is  determined  in  the  model  by  a 
parameterization  of  the  ageostrophic  mass  transport  in  the 
planetary  boundary  layer  (PEL) .  The  equation  relating  the 
horizontal  divergence  and  the  turbulent  Reynolds  shear 
stress  in  the  PBL  may  be  written  (  following  Haltiner,  1971) 
as: 


(2.12) 


where  T: 


d  are  the  horizontal  shear  stress  components 


EX  an 


in  the  x  and  y  directions  respectively  due  to  the  vertical 
wind  shear.  These  are  assumed  to  vanish  at  the  top  of  the 
PEL.  If  this  equation  is  integrated  from  the  bottom  to  the 
top  of  the  PBL,  the  resulting  expression  for  the 
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f rictionally -induced  vertical  velocity  component  CJ ^  at  the 
top  of  the  PBL  is: 

=  (!?■’-  <2-13' 

The  surface  stress  T  is  assumed,  in  this  model,  to 
be  proportional  to  the  square  of  the  wind  speed  \/0  at  a 
height  of  ten  meters  (anemometer  level)  according  to  the 
relationship : 


f  =  CA  V0  |l(|  (2-14) 

The  drag  coefficient  C £  depends  on  the  roughness  of  the 
underlying  surface,  while  P^  is  the  air  density  at  terrain 
height  Ht  .  As  with  the  pressure,  the  density  at  terrain 
height  is  given  the  value  corresponding  to  that  in  the  U.  S. 
Standard  Atmosphere.  Thus  the  frictional  vertical  velocity 
is  now  written: 


^  =3/f  klteCaT'.IV.l)-^  (ftQu0|\/J)]  (2.15) 


where  U0  and  V0  are  respectively  the  x-and  y-coraponents 
of  the  horizontal  wind. 


There  are  several  methods  by  which  V0  may  be 
estimated.  Some  quasi-geostrophic  models  use  the  value  of 
V/g0  at  a  level  near  terrain  height.  Danard  (1966)  for 
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example  used  the  850  mb  geostrophic  wind.  Some  models  use 
the  geostrophic  1000  mb  wind  forecast  by  the  model  while 
others  employ  a  wind  which  is  extrapolated  to  terrain  height 
from  the  flow  at  higher  levels  (e.g.  Cressman,  1963).  In 
the  model  used  in  the  current  study,  experiments  were 
conducted  using  the  850  mb  geostrophic  wind,  as  well  as  a 
geostrophic  wind  extrapolated  to  terrain  height  by  means  of 
the  Aitken  polynomial  which  is  described  in  Appendix  A. 


The  geostrophic  wind  V^0  may  be  further  modified  in 
the  calculation  of  the  vertical  velocity  CO^  .  Following 
Sawyer  (1959),  Greystone  (1962)  and  Haltiner  (1971),  a 
representative  value  for  X/0  may  be  considered  to  be 
smaller  by  a  fraction  r,  than  the  magnitude  of  the 
geostrophic  wind  \/<jo  at  the  top  of  the  planetary  boundary 
layer.  This  wind  may  also  be  turned  through  some  angle 
from  the  direction  of  \4j0  .  Making  these  assumptions  it 

may  be  shown  from  a  consideration  of  Figure  3  that: 


U„  =■  r  (u3o  Cos  oL  -  V%<3  S  IN  ) 

Vo  =  r  (y%o  *t-Ug0 S IN <*) 


(2.  16) 


where  and  V^0 are  respectively  the  x-and  y- components  of 

the  geostrophic  wind  Vg0 .  Substituting  into  equation 
(2.15),  the  complete  expression  for  the  friction ally- induced 
vertical  velocity  at  the  top  of  the  boundary  layer  becomes: 
CO  f  -  $/f  Cl  Y  (l^0C05  -f  IL^0  S  INo(  )  |  V<jJ  ] 


'4 


{ft0*  ^(V05^^51^ 


)|V,.lU 


(2.17) 


17 


Figure  3.*  The  relationship  between  the  geostrophic  wind 
vector  Vg0  at  the  top  of  the  planetary  boundary  layer, 
surface  stress  and  the  resultant  wind  vector  \4,  at 
anemometer  level  employed  in  the  calculation  of  the 
f rict ionally-induced  vertical  velocity  at  the  top  of  the 
PBL. 


the 
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Both  the  reduction  factor  r  and  the  turning  angle  of 
are  dependent  on  the  nature  of  the  underlying  surface. 

Sawyer  (1959)  suggests  that  may  vary  from  6°  over  the 

sea  to  30°  over  land  and  to  35°  over  mountainous  areas. 
Greystone  (1962),  in  his  two-level  model,  considered  values 
for  ^  of  3°  to  6°  over  the  sea  and  30°  elsewhere. 

The  value  of  r  used  in  Greystone® s  model  was  .35  over 
land  and  .85  over  the  sea,  while  the  drag  coefficient  had 
one  constant  value  over  land  and  another  value  over  the  sea. 
In  the  models  of  Cressman  (1960,1963),  Danard  (1966)  and 
many  others,  the  drag  coefficients  were  variable  in  space, 
with  a  maximum  value  over  mountain  ridges  and  a  smaller 
constant  value  over  the  sea  surfaces  and  flat  land.  The 
drag  coefficients  used  in  the  current  study  are  the 
empirical  values  calculated  by  Cressman  (1960),  which  are 
modified  by  the  inclusion  of  a  variable  r  and  as 

indicated  in  equation  (2.17).  The  justification  for  such  a 
modification  is  that  we  are  attempting  to  determine  the 
parameterization  of  the  lower  boundary  condition  which  best 
simulates  atmospheric  development  near  the  mountains.  A 
discussion  of  the  different  combinations  of  the  parameters 
used  at  the  lower  boundary  will  be  given  in  Chapter  V. 

While  is  assumed  to  originate  at  the  top  of  the 

PBL ,  in  most  filtered  equation  models  with  a  pressure 
coordinate  system,  the  thickness  of  this  layer  is  ignored  in 
comparison  with  the  thickness  of  the  next  layer  above. 


. 
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Thus,  in  these  models  is  considered  to  originate  at 

terrain  height.  If  the  thickness  of  the  PBL  were  taken  into 
account,  there  would  be  two  options  available:  the  PBL 
thickness  could  be  specified  as  constant  in  time  and 
variable  in  space,  or  a  predictive  equation  for  this  layer 
could  be  specified.  The  first  option  is  equivalent  to  a 
modification  of  the  terrain  height  field  which  does  not  seem 
justified  because  of  the  smoothed  topography  already  present 
in  the  model.  The  second  option  is  difficult  to  employ  in  a 
filtered  model,  although  is  is  more  amenable  to  inclusion  in 
a  primitive  equation  model  using  the  sigma  pressure 
coordinate  system. 

Thus  in  filtered  models,  both  and  are  assumed 

to  originate  at  terrain  height  Ht  (and  hence  )  .  An 

even  more  crude  assumption  often  used  for  convenience  is 
that  the  terrain-induced  vertical  velocity  originates  at 
some  lower  pressure  surface,  for  example  the  1000  millibar 
surface.  Such  an  assumption  allows  a  simpler  representation 
of  vertical  derivatives  near  the  lower  boundary.  A  full 
description  of  the  different  boundary  conditions  tested  in 
the  model  will  be  presented  in  Chapter  V. 
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CHAPTER  III 


DATA  ACQUISITION  AND  ANALYSIS 


3 • 1  H£D er  Air  Data 

For  two  of  the  three  cases  chosen  for  study  (January  5 
and  March  5f  1972) ,  the  upper  air  data  were  extracted  from 
the  Northern  Hemispheric  Data  Tabulations,  (NHDT) ,  available 
from  the  National  Climatic  Center  in  Asheville,  North 
Carolina,  U.S.A.  These  data  were  in  the  form  of  radiosonde 
reports  from  which  the  temperature,  wind,  and  geopotential 
height  at  the  four  primary  levels  were  coded  on  cards.  For 
about  85  upper  air  stations  lying  on  or  near  the  grid,  the 
above  data  were  extracted  for  both  the  initial  time  and  for 
the  verification  time,  12  hours  later.  In  the  third  case 
(April  27,  1974) ,  the  data  were  obtained  from  the  records  of 
the  Canadian  Meteorological  Center  (CMC)  at  Dorval,  Quebec. 
This  data  was  made  available  in  magnetic  tape  form,  and 
consisted  of  the  output  of  the  CMC  objective  analysis 
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program  on  the  381  km  hexagonal  grid  used  in  the  hemispheric 
numerical  weather  prediction  model  run  by  CMC. 

With  the  radiosonde  data  available  at  standard  upper 
air  stations  only,  a  large  area  of  the  grid  is  not  covered 
with  observations  in  the  first  two  cases.  These  data  sparse 
areas  over  the  Eastern  Pacific  Ocean  are  especially 
troublesome  when  attempting  objective  analysis,  because  the 
computer  program  described  in  Section  3.3  will  extrapolate 
unrealistic  flow  patterns  and  gradients  into  these  areas. 
Consequently,  bogus  data  must  be  generated.  A  sufficient 
density  of  bogus  data  points  were  selected  so  that  no  grid 
point  was  further  than  2  or  3  grid  intervals  from  any  bogus 
data  point.  Geopotential  heights,  calculated  geostrophic 
winds  and  temperatures  were  then  estimated  at  these  points 
from  the  CMC  objective  analysis  available  in  map  form. 

Since  these  maps  are  based  on  an  analysis  with  a  much  more 
extensive  data  base  than  that  available  in  this  study,  a 
consistent  and  reasonably  accurate  analysis  should  thereby 
be  possible  in  the  data  sparse  areas  of  this  study.  No 
testing  of  alternate  bogussing  techniques  was  undertaken. 

Since  no  300  mb  maps  were  available,  it  was  necessary 
to  extrapolate  the  300  mb  bogus  data  from  that  at  700  and 
500  mb.  For  this  purpose  a  simple  equation  was  obtained  by 
integration  of  the  hydrostatic  equation  from  700  and  500  mb 
to  the  300  mb  pressure  surface.  The  resulting  equation  is: 
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Z 3={ls+7.57{Js  -  r?)-0.6Z7]/ 0.3  9  S’  f3-1) 

where  Z7,Z  and  £^3  are  the  geopotential  heights  in 
meters  at  700,  500  and  300  mb  respectively,  and  T7  and 
Ts-  are  the  temperatures  in  degrees  Celsius  at  700  and  500 
mb.  This  equation  was  tested  at  radiosonde  stations  in  the 
vicinty  of  the  data  sparse  area,  and  a  representative 
correction  value  was  determined  based  on  the  mean  difference 
between  the  reported  and  extrapolated  300  mb  height  at  these 
stations.  This  corrected  extrapolation  equation  was  then 
used  to  generate  300  mb  heights  over  the  Eastern  Pacific 
portion  of  the  grid. 

3 • 2  Surface  Data 

The  terrain  heights  used  in  the  model  were  those 
extracted  by  Schram  (1974).  The  objective  technique  for 
obtaining  these  heights  at  grid  points  consisted  in 
determining  the  average  terrain  height  within  a  circle  of 
radius  100  km  around  each  grid  point  using  1:1,000,000  World 
Aeronautical  Charts.  In  the  current  model,  this  terrain 
height  field  is  converted  to  pressure  in  the  U.S.  Standard 
Atmosphere  (Hess,  1959)  by  making  use  of  the  formula 
Pt  =1  01  3 . 25[  1-  (C.  0065/28  8)  Ht  ]S,2S7.  A  similar  formula 
was  used  to  determine  the  atmospheric  density  at  terrain 
height.  Figure  4  depicts  the  terrain  pressure  field  used  in 
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Figure  4.  Pressure  at  terrain  height  calculated  using  the 
United  States  Standard  Atmosphere.  The  units  are  millibars. 
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Figure  5.  The  analysis  of  the  Cressman  drag  coefficients. 
The  units  are  10~3  (dimensionless). 
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the  model. 

The  drag  coefficients,  which  were  supplied  by 
Dr.  G.  Cressman  to  Hr.  G.  Schram,  have  been  interpolated  to 
the  200  km  grid  used  in  this  study.  The  drag  coefficients 
were  in  the  form  of  grid  point  data  as  discussed  by  Cressman 
(1960)  in  his  study  with  the  barotropic  model.  These 
coefficients  are  composed  of  two  components:  a  constant 
value  of  1.296x10~3  (dimensionless)  over  all  surfaces,  and 
an  additive  quantity  which  is  a  function  of  terrain 
roughness.  Thus  the  drag  coefficient  ranges  up  to  slightly 


UNDERLYING  SURFACE 
WATER 
LAND 

MOUNTAINS 


103  X  DRAG  COEFFICIENT 
Cd<1 .3 
1 .3<Cd<4. 0 
Cd>  4.0 


Table  1.  Surface  regimes  as  a  function  of  drag  coefficient. 


in  excess  of  8x10”3,  the  latter  over  very  rough  mountains. 
Figure  5  shows  the  drag  coefficient  over  the  area  of  this 
study. 

As  indicated  in  the  discussion  of  the  lower  boundary 
condition,  the  turning  angle  ,  o(  ,  and  the  reduction 
factor,  r,  are  variable  in  space  depending  on  the  nature  of 
the  underlying  surface.  In  this  study,  three  regimes  were 
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chosen,  in  each  of  which  o4  and  r  assume  constant  values. 
The  regimes  represent  mountainous  areas,  relatively  level 
land  surfaces  and  extensive  water  surfaces.  The  criterion 
for  distinguishing  among  these  areas  is  based  on  the 
magnitude  of  the  Cressman  drag  coefficients.  In  this  study, 
the  three  regimes  and  the  corresponding  range  of  drag 
coefficients  are  given  in  Table  1.  The  particular  values  of 

and  r  tested  are  discussed  in  Chapter  V. 

3-3  The  Objective  Analysis 

An  objective  analysis  program  similar  to  that  described 
by  Glahn  and  Hallenbaugh  (1969)  is  used  to  interpolate  the 
meteorological  fields  to  the  grid  points  for  machine 
calculation.  Studies  of  a  similar  type  of  objective 
analysis  procedure  by  Stuart  (1974)  suggest  that  such  a 
technique  compares  favorably  with  the  most  careful  hand 
analysis.  The  objective  method  is  essentially  one  of 
successive  approximation,  and  is  completely  described  by 
Glahn  and  Hallenbaugh  as  well  as  by  Bergthorssen  and  Doos 
(1955)  and  by  Cressman  (1959).  It  is  the  technique 
frequently  applied  in  the  analysis  of  spatial  fields  for 
numerical  weather  prediction  and  analysis. 

Since  small  scale  irregularities  in  the  analyzed  fields 
may  amplify  during  time  integration  of  the  numerical  model, 
it  is  desirable  to  remove  such  two- grid-inter val 
perturbations  from  the  initial  fields.  All  height  fields 
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were  therefore  smoothed  in  the  first  two  cases  using  the 
techniques  described  in  Appendix  B.  The  resultant  fields 
appeared  reasonably  accurate  over  the  continent.  Rhen 
compared  with  the  reported  radiosonde  data,  the  analyzed 
fields  had  a  mean  absolute  deviation  over  the  entire  grid  of 
approximately  20  gpm  at  all  four  levels.  The  accuracy  over 
the  data  sparse  areas  is  of  course  difficult  to  assess  and 
hence  only  a  qualitative  judgement  could  be  made.  However, 
with  the  use  of  a  sufficient  density  of  bogus  data  points, 
as  suggested  above,  the  analysis  was  consistent  with  that  of 
the  CMC.  At  300  mb  no  such  comparison  could  be  made  and,  as 
described  in  Chapter  V,  some  difficulty  with  prediction  at 
this  level  may  have  been  partly  due  to  an  incorrect 
initialization  of  the  geopotential  field. 

In  the  third  case,  the  upper  air  geopotential  height 
data  were  extracted  from  a  magnetic  tape  supplied  by  CMC.  A 
simple  program  was  written  to  extract  the  data  from  a 
section  of  the  CMC  grid  which  contained  the  smaller  grid 
used  in  this  study.  Since  the  CMC  grid  spacing  is  381  km, 
it  was  necessary  to  perform  an  interpolation  to  the  200  km 
grid  of  this  study  utilizing  the  objective  analysis  program. 
No  smoothing  or  bogussing  of  the  resultant  geopotential 
fields  was  done,  as  the  original  CMC  analysis  was  heavily 
smoothed. 

Since  the  surface  pressure  (mean  sea  level  pressure) 
field  is  not  forecast  directly,  the  actual  surface  pressure 
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data  does  not  enter  into  the  mode]..  However,  since  it  is 
desirable  to  obtain  a  forecast  surface  pressure  map  for 
comparison  with  the  actual  surface  analysis,  such  a  field 
was  obtained,  using  the  Aitken  polynomial  to  extrapolate  the 
surface  pressure  from  the  forecast  geopotential  fields 
aloft.  In  order  to  assess  the  errors  which  could  be 
expected  from  such  a  procedure,  the  surface  pressure  field 
at  the  initial  time  in  the  first  case  was  objectively 
analyzed  using  a  program  similar  to  that  already  described. 

A  comparison  of  this  actual  mean  sea  level  pressure  field 
and  the  extrapolated  pressure  field  appears  in  Chapter  V. 
Generally  there  was  good  agreement  between  the  extrapolated 
and  the  actual  fields  when  the  broad  scale  features  were 
compared. 

3 • 4  The  Static  Stability  Parameter 

The  static  stability  parameter  o'  used  in  the  model 
is  a  function  of  pressure  only,  and  is  computed  using  the 
temperatures  available  in  the  radiosonde  ascents.  A  simple 
second-order  centered-dif f erence  approximation  to  the  static 
stability  equation  O'  =  -  ~  was  formulated  and,  for 

each  of  the  three  intermediate  levels.  O'  was  calculated  at 
each  radiosonde  station  in  the  forecast  area.  A  horizontal 
average  O'  was  then  calculated  at  each  level  and  used  in 
the  solution  of  the  omega  equation.  Since  a  large  station- 
to-station  variability  was  noted  in  O'  ,  the  standard 
deviation  at  each  level  was  also  calculated.  The  numerical 
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values  of  the  mean  stability  parameter  and  the  standard 
deviation  are  given  in  Table  2.  In  general,  the  largest 
absolute  deviations  from  the  mean  occur  at  high  latitudes  at 
the  400  mb  level,  since  the  tropopause  is  approaching  this 
level  at  these  latitudes.  At  low  levels,  there  are  a  few 
stations  indicating  an  unstable  stratification  (O'  <0)  .  Thus 
unfortunately,  the  mean  values  used  are  not  representative 
of  any  particular  region  of  the  grid. 

In  the  current  model,  some  testing  of  different  values 
of  the  static  stability  parameter  was  performed.  A  mean 
stability  was  calculated  for  only  those  radiosonde  stations 
lying  in  the  section  of  the  grid  between  45°  and  65°  north 
and  to  the  west  of  the  Manitoba- Saskatchewan  border.  This 
mean  was  found  to  differ  from  the  overall  grid  average  by  up 
to  twenty-five  percent.  However,  when  these  stabilities 
were  employed  in  the  solution  of  the  omega  equation,  only 
very  small  percentage  changes  (generally  less  than  five 
percent)  in  the  vertical  velocity  and  height  tendency  fields 
were  noted.  Comparisons  of  the  vertical  velocity  and 
forecast  height  fields  as  calculated  using  variable  and 
constant  static  stabilities  by  Haltiner  et  al  (1963)  and  by 
Cressman  (1963)  have  indicated  that  only  small  percentage 
changes  occured. 


Since  it  appears  that  the  solution  to  the  omega 
equation  is  relatively  insensitive  to  the  manner  in  which 
the  static  stability  is  included  in  numerical  models,  it 
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775  mb 

CASE  MEAN  SD 

1  2.  55  2.  84 

2  2.21  2.38 

3  1.79  1.92 


600  mb 
MEAN  SD 
2.38  2.52 
2.32  2.41 
2.14  2.17 


400  mb 
MEAN  SD 
4.93  5.31 
4.60  5.11 
3.79  4.68 


Table  2.  Mean  and  standard  deviation  (SD)  of  the  static 
stability  over  the  complete  grid  (10~4  gra~4  cm4  s2)  at  three 
levels  in  the  model. 


would  seem  justifiable  to  use  a  constant  value  of  the 
stability  at  each  level.  In  fact  many  models  have  employed 
a  stability  based  on  mean  yearly  values  or  on  those 
calculated  from  a  standard  atmosphere. 
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CHAPTER  IV 


NUMERICAL  PROCEDURES 


4 . 1  Horizontal  Finite-Dif fere nee  Operators 

a)  First  Derivatives 

The  standard  centered-dif ference  formulation  is  used 
for  the  horizontal  derivatives  appearing  in  equation  (2.17), 
which  is  the  expression  for  the  vertical  velocity  at  the 
lower  boundary.  Thus,  with  reference  to  Figure  6  which 
illustrates  a  nine  point  section  of  the  interior  of  the 
grid,  the  x-and  y-der ivati ves  of  the  field  A  at  the  point 
i,j  are  represented  by  the  second-order  formulae: 


(4.1) 


■ 


32 


Y 


i+l»j  “I 


i+l’j 


i+lj+l 


•  •  *  •  • 

l.J-1  !.J 


i.j  +  1 
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i  l»j  +  1 


Figure  6.  The  horizontal  grid  point  structure  used  in  the 
finite  differencing.  The  Y-axis  corresponds  to  true  north 
at  longitude  115°  W. 


33 


i+ 1  j 


(4.2) 


At  the  boundaries  of  the  grid,  the  horizontal 
derivatives  are  approximated  by  one-sided  derivatives.  For 
example  if  point  i+1,j  lay  along  the  northern  edge  of  the 
grid,  then  the  x-derivative  would  be  calculated  in  the  usual 
manner,  whereas  the  y-der i vative  would  be  approximated  by 
the  first-order  formula: 


(4.3) 


A  similar  technique  is  used  to  calculate  the  x-derivative  at 
the  eastern  and  western  boundaries  and  the  y-derivative  at 
the  southern  boundary. 

As  indicated  in  Section  2.5,  the  horizontal  wind  used 
in  the  calculation  of  the  terrain-induced  vertical  velocity 
is  a  function  of  the  geostrophic  wind.  In  the  model,  this 
may  be  either  the  850  mb  geostrophic  wind  or  a  geostrophic 
wind  extrapolated  to  terrain  height.  Since  the  850  mb  level 
is  a  forecast  level,  the  calculation  of  the  wind  in  the 
first  instance  is  readily  performed  by  taking  the  horizontal 
derivatives  of  the  850  mb  height  field. 


In  order  to  obtain  extrapolated  winds,  the  information 
at  the  four  primary  levels  in  the  vertical  is  extrapolated 
to  terrain  height  by  making  use  of  the  Aitken  polynomial. 
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While  the  geostrophic  winds  calculated  at  each  of  these 
levels  above  a  grid  point  i,j  could  be  extrapolated  to 
terrain  height  directly,  it  was  found  that  rather 
unacceptable  winds  would  be  extrapolated  if  large  vertical 
wind  shears  occured.  This  problem,  inherent  in  the  use  of 
the  Aitken  polynomial,  is  greatly  diminished  if  the 
geopotential  height  is  extrapolated,  and  then  the 
geostrophic  wind  calculated.  Thus  with  respect  to  point  i,j 
of  Figure  6,  the  geopotential  height  is  extrapolated  to  an 
isobaric  surface  which  intersects  the  terrain  at  pressure 

.  The  height  at  each  of  the  four  points  i,j±1  and 
i±1,j  on  this  isobaric  surface  is  then  extrapolated  from  the 
geopotential  height  at  the  primary  levels  above  each  of 
these  points.  Thus  four  extrapolations  per  grid  point  are 
required  in  order  to  calculate  both  components  of  the 
geostrophic  wind  vector  at  point  i,j. 

The  procedure  is  repeated  for  each  of  the  19x21 
interior  grid  points  and  allows  a  second-order  calculation 
of  the  geostrophic  wind  to  be  made  at  terrain  height.  Along 
the  edges  of  the  grid,  extrapolation  is  made  in  a  similar 
manner,  except  that  one-sided  first-order  derivatives  are 
calculated  for  the  wind  components  parallel  to  the  boundary 
of  the  grid.  In  cases  where  the  terrain  height  rises  above 
the  850  mb  surface,  an  interpolation  is  performed  making  use 
of  the  geopotential  heights  of  the  four  primary  levels. 

b)  The  Two-dimensional  Laplacian 


' 
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With  reference  to  Figure  6,  the  standard  finite- 
difference  five-point  Laplacian  operator  for  variable  A  may 
be  expressed  as: 


(7‘A)  =(A,.litAt.11+All^AtJ.1-1f-Als)/a-r^5  <4-4> 


d.-l 


This  operator  is  accurate  to  order  d2(-  . 


c)  The  Jacobian 


With  the  Jacobian  of  two  variables  A  and  B  defined  as 
J (A , B)  = JA ^3- •  the  finite-difference  equivalent  may  be 
written  in  its  five-point  form  as: 


J  (ABXj  =  [(Aij+1  ACj_,)(  Bu.j  B^.j) 


■+  + 


(4.5) 


-(Btj4rBij-.)(AUirA.l.li)]/fdL*j=J‘ 

where  J”  is  the  finite- difference  representation  for  X  . 


In  prediction  models,  non-linear  computational 
instability  may  arise  if  the  representation  of  the  advective 
(or  Jacobian)  terms  in  the  predictive  equation  is  not 
performed  so  as  to  conserve  the  mean  square  vorticity  and 
the  energy  of  the  system.  Arakawa  (1966)  has  shown  that  a 
finite-difference  Jacobian  of  the  form: 


J  -  (X+++  Jx%  J+*)/3 


will  conserve  both  energy  and  mean  square  vorticity.  The 
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various  Jacobian  operators  in  this  equation  are  defined  as 
follows: 


-r/A  p\  <2)  /A  8  A  AG 

-J  V'v  D  J  —  x  <^x 

=  J++ 

a 

A  /  a  A  G  ^  _  A  /  a  A  G  'j 

=  J+x 

lr  (4.7) 

-  ^  / g ^  \ _ A  (B  —  ) 

^Jx+ 

c 

The  f inite-di f ference  formulations  analagous  to  equation 

(4.5)  for  each  of  b  and  c  in  equation  (4.7)  may  then  be 
combined  into  the  nine-point  Arakawa  Jacobian  for  use  in  the 
vorticity  equation. 

4.2  The  Vertical  Derivatives 

Since  this  model  has  unequally  spaced  information 
levels  in  the  vertical,  it  is  not  possible  to  use  a  simple 
centered-dif ference  formulation  for  the  first  and  second 
derivatives  of  co  which  appear  in  equations  (2. 3)  and 

(2.5) .  Instead  a  weighted  difference  scheme  is  used. 
Haltiner  et  al  (1963)  calculated  these  derivatives  based  on 
the  assumption  of  a  quadratic  variation  of  CO  with  height 
over  three  successive  levels.  In  this  study,  a  similar 
overlapping  parabolic  profile  is  assumed  with  the  finite- 
difference  formulation  based  on  the  second  order  Aitken 
polynomial.  The  finite-difference  approximations  for  the 
first  and  second  derivatives  of  the  Aitken  polynomial  are 
discussed  in  Appendix  A. 
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In  several  portions  of  the  grid,  the  terrain  rises 
above  the  850  mb  surface  to  just  below  the  700  mb  surface. 

In  these  areas,  special  consideration  has  to  be  given  to  the 
evaluation  of  the  vertical  derivatives  of  the  CO  field. 
While  it  may  not  be  physically  meaningful  to  evaluate  the 
forcing  functions  in  the  vorticity  equation  at  the  850  mb 
level  if  that  surface  lies  below  the  terrain  surface, 
changes  in  the  vorticity  and  flow  patterns  at  850  mb  still 
occur  due  in  part  to  the  influence  of  the  topography.  Thus 
it  is  desirable  to  parameterize  these  changes  even  if  the 
terrain  in  places  lies  slightly  above  the  850  mb  surface. 

In  portions  of  the  grid  where  the  terrain  rises  above  850 
mb,  the  divergence  at  850  mb  was  extrapolated  from  that 
calculated  at  the  terrain  height.  The  820  mb  level  was 
chosen  as  the  cutoff  level  beyond  which  no  extrapolation  was 
attempted.  This  surface  lies  just  above  the  highest 
mountain  ridges  in  the  central  and  northern  portions  of  the 
grid.  Where  the  terrain  rises  above  the  820  mb  level,  the 
850  mb  divergence  was  set  to  zero.  Thus  at  850  mb,  the 
horizontal  divergence  was  calculated  over  all  of  the  region 
of  interest  even  though  the  terrain  rises  above  the  850  mb 
surface  at  several  points.  Since  the  terrain  does  not  rise 
above  the  700  mb  level,  no  problems  arose  in  the  evaluation 
of  the  divergence  at  that  level. 

In  the  omega  equation,  the  second  derivative  of 
the  CO  field  with  respect  to  pressure  is  evaluated  at  the 
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775  mb  level.  In  regions  of  high  terrain,  this  term  was 

evaluated  in  a  manner  consistent  with  the  calculation  of  the 

first  derivative  of  6J  .  Since  by  the  continuity 

equation,  the  second  derivative  of  GJ  with  respect  to  p  is 

equivalent  to  the  rate  of  change  of  horizontal  divergence 

with  height,  the  divergence  both  above  and  below  775  mb 

^U> 

enters  into  the  calculation  of  •  In  regions  where  the 

terrain  lay  above  820  mb,  the  evaluation  of  the  rate  of 
change  of  divergence  with  height  involved  only  the 
divergence  above  775  mb,  the  divergence  below  this  level 
being  set  equal  to  zero. 

4 . 3  The  Time  Step 

The  vorticity  equation  used  in  this  study  is 
essentially  a  non-linear  advection  equation  with  a  sink  or 

r  ^  6J 

source  of  vorticity,  T  §"p  .  Omitting  this  sink  term,  a 
linearized  one-dimensional  form  of  the  vorticity  equation 
is: 

=  -cfx  (<r+f)  0-8) 

where  C  is  the  advection  speed  of  the  vorticity  C  +  f  in  ihe 
x- direction.  If  the  advection  term  on  the  right- hand-side 
of  this  equation  is  approximated  by  a  centered- difference 
formulation,  and  a  centered  time  difference  is  employed, 
then  the  Courant-Friedrichs-Lewy  (CFL)  condition  for  linear 
computational  stability  of  the  numerical  solution  of 


. 

■ 
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equation  (4.8)  is: 


At  £ 


A  X 

C 


(4.9) 


on  a  grid  with  spacing  Ax  employing  a  time  step  of  length 
At.  The  two  dimensional  analogue  of  equation  (4.9)  can  be 
shown  to  be: 


At 


AX 

•Iz'c 


(4.10) 


In  filtered  models,  in  which  rapidly  propagating  gravity  or 
sound  waves  cannot  occur,  C  will  be  the  maximum  wind  speed 
in  the  atmosphere.  In  the  cases  chosen,  this  maximum  speed 
occurs  at  the  300  mb  level,  and  with  the  grid  spacing  of 
approximately  200  km  used  in  this  model,  the  maximum 
allowable  time  step  is  in  the  range  .5  to  .75  hours 
depending  on  the  300  mb  wind  speed. 

While  the  correct  choice  of  the  time  increment  is  a 
sufficient  condition  to  ensure  linear  computational 
stability,  a  non-linear  instability  may  be  associated  with 
the  finite-difference  representation  of  the  advection  term. 
As  indicated  in  Section  4.1,  the  Arakawa  Jacobian,  since  it 
conserves  both  mean  square  vorticity  and  mean  kinetic  energy 
should  eliminate  this  nonlinear  instability.  However,  as 
Haltiner  (1971)  points  out,  these  conservative  properties 
have  only  been  proven  for  the  continuous  time  derivative 
case  and  not  for  the  finite-differencing  time  step. 


:r\  *« 
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The  standard  centered-time-dif ference  technique  used  to 
forecast  a  new  field  X  is: 


(4.11) 


where  n  is  an  integer  indicating  the  number  of  time  steps 


from  the  beginning  of  the  forecast,  and 


tendency  field,  which  in  the  current  model  is  calculated  by 
numerical  representation  of  the  right-hand-side  of  the 
vorticity  equation.  The  procedure  is  started  with  a  forward 
time  step  according  to: 


(4.12) 


As  shown  by  Haltiner  (1971),  when  this  f orward- then- 
centered-dif ference  technique  is  employed,  neither  the  the 
mean  square  vorticity  nor  mean  energy  of  the  flow  are 
conserved  even  with  the  Arakawa  formulation  of  the  Jacobian. 
The  result  is  a  non-linear  instability,  which  is  initially 
manifested  as  a  small  amplitude  temporal  oscillation  of  the 
numerical  solution. 

As  indicated  by  Gerrity  (1972),  this  oscillation  or 
separation  of  the  solutions  between  the  even  and  odd  time 
steps  is  generally  not  a  serious  problem  in  most  numerical 
models  provided  that  the  variation  of  the  vorticity  field 
between  successive  time  steps  is  small.  However,  in  a 


1 1 . 
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baroclinic  model  such  as  that  of  the  current  study,  the 
sources  of  vorticity  which  arise  from  the  flow  of  the 
atmosphere  over  steep  terrain  may  show  large  spatial  and 
temporal  variations.  The  resulting  separation  of  the 
numerical  solutions  with  time  then  becomes  very  large,  and 
at  some  point  in  time  the  oscillation  reaches  such  a 
magnitude  that  the  predicted  fields  are  severely  distorted 
and  further  prediction  is  meaningless.  This  is  the 
manifestation  of  non-linear  instability. 

This  difficulty  has  been  recognized  for  some  time,  and 
many  procedures  have  been  suggested  to  overcome  this  so- 
called  ’weak  instability’.  Gates  (1960)  used  a  technique  in 
which  the  forecast  was  begun  with  a  short  forward  time  step 
followed  by  a  sequence  of  centered  time  steps  of  increasing 
length  until  the  final  desired  length  of  time  step  was 
achieved.  Such  a  procedure  reduces  the  inital  amplitude  of 
the  oscillations  but  does  not  filter  them  out  completely. 
Thus  the  oscillation  still  increases  in  time  and  the 
procedure  has  to  be  re-initiated  after  several  time  steps 
following  some  temporal  smoothing. 

An  investigation  of  several  time  stepping  techniques 
has  been  undertaken  by  many  authors,  especially  with 
reference  to  the  primitive  equation  mode  s  in  which  some 
temporal  smoothing  is  generally  found  to  be  necessary.  In  a 
paper  by  Lilly  (1965) ,  several  finite-difference  operators 
were  tested  in  detail  for  a  system  of  advection  equations. 


iv  m  Ml 
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The  numerical  and  analytical  solutions  to  the  simplified 
system  of  equations  were  compared,  and  it  was  shown  that  the 
standard  cen  tered-dif ference  formulation  in  time  exhibited 
extreme  instability  resulting  in  time  splitting  of  the 
numerical  solution  when  certain  combinations  of  initial 
conditions  were  specified.  Several  other  time  stepping 
techniques  of  varying  complexity  and  computational 
efficiency  were  also  tested,  and  all  appeared  to  have 
improved  stability.  Lilly  indicated  that  the  Adams- 
Bashforth  method  seemed  to  be  one  of  the  most  efficient  in 
terms  of  computation  time  required  per  time-step,  and  that 
it  exhibited  close  agreement  with  the  analytical  solution  to 
the  simplified  equations.  The  difference  formulation  for 
this  technique  is: 


(4.13) 


where  again  n  indicates  the  time  level,  X  is  the  predicted 


field  and  T  the  tendency  field. 


Further  tests  on  the  accuracy  of  various  time  stepping 
techniques  were  undertaken  by  Molenkamp  (1968)  again  using  a 
simple  advection  model.  A  comparison  of  several  time 
stepping  procedures  suggested  that  the  Adams- Bashforth 
technique  was  of  comparable  accuracy  to  the  centered- 
difference  formulation  when  both  were  used  to  predict  the 
advection  of  simple  patterns.  As  well,  the  Adams-Bashf orth 
technique  showed  a  smaller  tendency  towards  damping  of  the 
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aavected  fields  than  did  the  centered~d.if  fere  nee 
formulation. 

In  the  present  study,  forecasts  were  first  computed 
with  the  f irst-f orward- then-centered-dif f erence  technique, 
but  it  became  evident  that  a  tendency  toward  separation  of 
the  solutions  at  alternate  time  steps  increased  with  certain 
lower  boundary  conditions.  In  Figure  7,  the  differences 
forecast  from  the  initial  value  of  the  850  and  500  mb 
geopotential  height  at  grid  point  (12,14)  are  plotted  as  a 
function  of  time.  The  grid  point  chosen  is  the  closest  one 
to  Edmonton,  Alberta  as  shown  in  Figure  1.  The  oscillatory 
nature  of  the  solution  is  most  apparent  at  850  mb  where  the 
lower  boundary  condition  is  most  influential  and  the 
amplitude  of  the  oscillation  rapidly  increases  with  time. 

By  about  nine  hours  (or  18  time  steps) ,  the  oscillations 
have  commenced  at  higher  levels  and  the  forecast  has  long 
since  become  meaningless.  When  the  Adams-Bashf orth  time 
step  is  used  and  all  other  initial  conditions  are  unchanged, 
the  forecast  shows  no  tendency  towards  oscillation.  Thus 
the  Adams-Bashforth  time  step  yields  a  computationally 
stable  solution  to  the  time-dependent  equations  for  this 
case  and  no  further  problems  attributable  to  the  time 
stepping  procedure  were  encountered. 

4.4  Solution  of  the  Differential  Equation s 


The  right  hand  sides  of  the  vorticity  and  omega 
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Figure  7.  Forecast  height  differences  from  the  initial 
condition  as  a  function  of  time  at  850  mb  and  at  500  mb:  (a) 
forward  then  centered  time  step,  (b)  a  single  forward 
followed  by  Adams-Bashf or th  time  step.  Both  were  calculated 
at  grid  point  12,14  which  is  indicated  in  Figure  1. 
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equations  consist  of  forcing  functions  composed  of 
differential  operators  acting  on  other  differential 
operators;  for  example,  the  Jacobian  operating  on  the 
Laplacian  of  the  height  field.  Because  of  the  manner  in 
which  the  operators  are  formulated  in  the  finite-difference 
scheme,  these  forcing  functions  can  be  evaluated  only  over  a 
portion  of  the  grid  which  excludes  a  strip  two  grid  points 
wide  all  around  the  lateral  boundary.  The  solution  to  the 
omega  equation  then  consists  of  finding  a  simultaneous 
solution  to  a  set  of  969  finite-difference  equations  -  one 
equation  for  each  of  17x19x3  points  in  the  three  dimensional 
grid.  The  vorticity  equation  is  solved  at  each  of  four 
levels  independently,  and  therefore  a  solution  to  four  sets 
of  323  (17x19)  finite-difference  equations  is  required. 

The  solutions  to  the  finite-difference  forms  of  the 
vorticity  and  omega  equations  are  achieved  by  means  of  an 
iterative  process  known  as  Liebmann  over-relaxation.  This 
procedure  is  frequently  used  in  numerical  weather 
prediction,  and  the  technique  is  described  in  texts  by 
Thompson  (1961)  and  Haltiner  (1971).  A  brief  description  of 
the  technique,  and  the  parameters  used  to  achieve  the 
numerical  solution  follows. 

The  vorticity  and  omega  equations  may  be  written  in  the 
form  of  a  Helmholtz  equation  as  follows: 


VZT  +  F(T)  =  G(x,,jP>t) 


(4.14) 


. 
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where  T=T(x,y,p,t)  is  the  CO  or  the  height  tendency 
field  and  the  forcing  function  is  G.  F  is  a  function  which 
involves  the  second  derivative  of  CO  in  the  omega 
equation,  while  F  is  zero  in  the  vorticity  equation.  When 
equation  (4,14)  is  represented  in  finite-difference  form, 
the  following  relation  is  obtained: 


K  FUl 


(4.15) 


<;jK  is  the  m*th  approximation  to  the  field  T  at 


point  (i, j , k)  and  is  the  five-point  Laplacian  operator 

defined  in  equation  (4,4) .  The  functions  K  and  K'  involve 
the  grid  spacing  and  the  static  stability,  and  R  is  the 
residual  which  is  to  be  made  as  small  as  possible,  A  new 
estimate  for  the  T  field  may  be  derived  from  the  previous 
estimate  using  the  equation: 


—r-~  tA  a-  I  .  m  o  r/1 

Tc jK  =  T£jk  +&  Kcj 


(4,16) 


In  this  equation  /2  is  the  relaxation  parameter,  which 
should  be  chosen  so  as  to  maximize  the  rate  of  convergence 
of  the  iteration  technique.  Through  experimentation  it  was 
found  that  in  the  solution  of  the  vorticity  equation,  a 
value  of  0,45  for  /3  allowed  a  convergent  solution  to  be 
achieved  in  25  to  35  iterations  at  each  time  step.  The  same 
value  of  /3  was  used  in  all  three  cases  tested. 
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In  the  solution  of  the  omega  equation,  it  was  found 
necessary  to  use  a  value  of  /3  in  Case  3  different  from 
that  used  in  Cases  1  and  2.  A  numerical  solution  to  the 
omega  equation  was  achieved  in  15  to  25  iterations  if  a 
value  of  =0.285  was  used  in  Cases  1  and  2  while  in  Case 

3,  a  value  of  /3  =0.245  was  the  optimum  value  used  to 
achieve  rapid  solution.  The  rate  of  convergence  in  the 
solution  to  the  omega  equation  was  very  sensitive  to  the 
relaxation  parameter.  An  increase  of  only  8%  in 
resulted  in  a  doubling  of  the  number  of  iterations  required 
to  solve  the  omega  equation  in  Case  3.  Both  the  40  and  the 
height  tendency  fields  were  calculated  to  an  accuracy  such 
that  the  largest  residual  was  about  1%  of  the  synoptic  scale 
value  to  be  expected  from  a  scale  analysis  of  the 
vorticities  and  omega  equations.  This  corresponds  to  an 
error  of  about  1  C~ 2  microbars  s~ 1  for  the  <50  field  and 
10~5  gpm  s”1  for  the  height  tendency  field. 

The  numerical  solution  described  above  is  essentially 
the  solution  to  a  boundary  value  problem  at  each  time  step. 
In  general,  the  values  of  the  function  T  in  equation  (4.14) 
are  unknown  along  the  lateral  and  vertical  boundary  and  must 
be  specified  as  a  function  of  time.  In  the  omega  equation, 
the  upper  and  lower  boundary  conditions  .aye  already  been 
described.  Along  the  lateral  boundaries,  CJ  is  set  to 
zero  at  each  grid  point  in  the  outermost  two  rows  of  points 
around  the  grid. 
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In  the  vorticity  equation,  the  lateral  boundary 

condition  for  the  height  tendency  field  may  be  specified  in 

a  number  of  ways.  As  in  the  omega  equation,  the 

specification  is  made  in  the  outermost  two  rows  of  points. 

In  many  filtered  equation  models,  a  clamped  lateral  boundary 

^  2 

condition  or  -c-j-  =  0  is  assumed.  This  condition  is 
reasonable  if  the  area  of  solution  for  the  vorticity 
equation  covers  a  hemisphere  so  that  the  lateral  boundary 
extends  into  subtropical  regions.  Synoptic  scale  changes  in 
the  flow  patterns,  and  hence  in  the  geopotential  field,  vary 
very  slowly  with  time  in  these  regions  compared  to  the 
changes  observed  in  middle  to  high  latitudes.  Thus  for 
integration  periods  of  from  one  to  three  days  over  a 
hemisphere,  the  clamped  boundary  condition  does  not  result 
in  serious  errors. 


If  the  integration  area  is  small  however,  the 
assumption  of  a  clamped  lateral  boundary  condition  will 
result  in  more  serious  error,  especially  if  the  integration 
period  is  more  than  a  few  hours.  In  order  to  minimize  this 
problem,  the  grid  is  made  large  enough  that  the  region  of 
interest  at  the  center  of  the  grid  is  far  removed  from  the 
lateral  boundaries.  For  example  in  the  models  of  Petterssen 
et  al  (1962),  Pedersen  (1  962),  Greystone  (1962)  and  Danard 
(1966),  the  lateral  boundaries  are  approximately  4000  to 
5000  km  from  the  interior  region  of  interest,  with  the 
result  that  the  clamped  boundary  condition  has  only  a  small 
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effect  on  the  behavior  of  synoptic  scale  motions  at  the 
center  of  the  grid.  In  the  current  model  however,  the 
boundaries  are  located  in  regions  where  the  atmosphere  is 
undergoing  rapid  circulation  changes  and  the  clamped 
boundary  condition  would  be  erroneous. 

Many  regional  models  have  been  formulated  in  which  a 
small  area  grid  of  short  grid  length  is  meshed  into  a 
portion  of  a  large  scale  grid.  The  lateral  boundary 
condition  on  the  meshed  grid  is  specified  by  the  interior 
solution  of  the  large  scale  model,  which  may  itself  cover  a 
hemisphere.  In  fact,  several  meshings  have  been  employed 
going  down  to  a  grid  size  as  small  as  1/8  of  that  of  the 
hemispheric  grid.  The  filtered  equation  models  of  Hill 
(1968)  and  Shapiro  and  0s  Brien  (1970)  are  two  examples  in 
which  such  a  procedure  has  been  employed. 

In  the  current  study,  since  it  was  not  feasible  to 
develop  a  hemispheric  four-level  model  to  provide  the 
lateral  boundary  condition,  it  had  to  be  specified  in  a 
different  way.  In  an  early  baroclinic  model  run  by  Thompson 
and  Gates  (1956),  the  lateral  boundary  condition  for  the 
tendency  field  was  specified  to  be  the  observed  24  hour 
height  tendency  averaged  over  a  period  centered  at  the  time 
of  initialization  of  the  forecast.  In  the  current  study,  a 
height  tendency  field  was  calculated  based  on  the  height 
difference  between  the  initial  height  field  and  that 
observed  12  hours  after  the  initial  time.  This  tendency  was 
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used  as  a  lateral  boundary  condition  for  the  solution  of  the 
vorticity  equation  at  each  time  step.  The  lateral  boundary 
condition  for  the  tendency  field  was  therefore  constant  with 
time. 

With  such  a  specification  of  the  lateral  boundary 
condition,  the  model  is  a  simulation  model  rather  than  a 
true  forecast  model.  However,  the  solution  in  the  interior 
of  the  grid  should  now  be  able  to  reflect  variations  in  the 
lower  boundary  condition  which  is  the  purpose  of  this  study. 
Information  gained  on  the  importance  of  the  lower  boundary 
condition  could  then  be  included  in  a  meshed  prediction 
model. 

4 • 5  The  Forecast 

With  the  initial  and  boundary  conditions  specified  in 
the  manner  described  above,  a  forecast  is  achieved  following 
a  straightforward  marching  procedure.  At  each  time  step, 
the  omega  equation  is  solved  subject  to  the  lateral  and 
vertical  boundary  conditions  indicated  in  previous  sections. 
At  this  point  the  vorticity  field  is  explicitly  smoothed 
using  the  smoothing  operator  described  in  Appendix  B.  This 
smoothing  should  help  to  control  linear  computational 
instability.  The  vorticity  and  omega  fields  are  then  used 
in  the  vorticity  equation  which  is  solved  subject  to  its 
lateral  boundary  conditions.  The  above  cycle  (solution  of 
the  omega  then  vorticity  equation)  is  repeated,  after 
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generation  of  a  new  geopotential  height  field  at  each  time 
step,  until  the  desired  length  of  forecast  is  achieved.  The 
vorticity,  height  and  CO  fields  at  any  time  step  can  be 
stored  or  printed  out  for  later  analysis. 

All  programs  were  written  in  the  Fortran  language  and 
compiled  using  the  Fortran  H  compiler.  The  compilation  and 
running  were  done  on  the  University  of  Alberta  IBM  360/67 
computer  and  required  about  1  minute  of  processing  time  per 
hour  of  forecast.  Because  of  the  large  number  of  logical 
tests  written  into  the  program,  the  program  is  rather 
inefficient,  and  it  is  expected  that  a  large  saving  in 
processor  time  could  be  achieved  were  the  programs  rewritten 


in  a  streamlined  form. 
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CHAPTER  V 


RESULTS  AND  DISCUSSION 


5 - 1  Preliminary  Comments 

Although  three  cases  were  chosen  to  demonstrate  the 
influence  of  the  lower  boundary  condition  in  the  numerical 
prediction,  most  of  the  experimentation  was  conducted  on  the 
January  5,  1972  case,  hereafter  referred  to  as  Case  1.  This 
was  a  case  of  major  cyclogenesis  in  the  lee  of  the  Canadian 
Rocky  Mountains,  most  of  the  activity  occuring  in  the  12 
hour  period  following  the  initial  time  of  the  forecast.  In 
Case  2  (Mar.  5,  1972)  and  Case  3  (Apr.  27,  1974),  actual  lee 
cyclogenesis  did  not  occur  but  it  appeared  that  the  Rocky 
Mountains  had  an  effect  on  the  movement  and  intensity  of  the 
existing  synoptic  features.  The  synoptic  situations  for 
Cases  1  to  3  are  described  fully  below. 


In  the  following  discussions,  most  of  the  emphasis  will 
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be  on  the  850  mb  geopotential  field  although  a  discussion  of 
the  surface  pressure  fields  is  presented  in  the  Case  1 
forecasts.  In  Cases  2  and  3,  since  the  surface  pressure 
patterns  were  very  similar  to  the  850  mb  patterns,  little 
additional  insight  into  the  validity  of  different  lower 
boundary  conditions  could  be  gained  from  examination  of  the 
surface  pressure  forecasts. 

5. 2  Testing  of  the  Numerical  Procedures 

In  the  initial  testing  on  Case  1 ,  numerous  forecasts  of 
length  1  to  3  hours  were  run  to  test  the  numerical 
procedures  employed  in  the  model.  The  response  of  the  model 
to  the  reduction  factor  and  turning  angle  used  to  obtain  the 
boundary  layer  wind  was  tested.  To  evaluate  the  influence 
of  different  combinations  of  and  r,  three  hour  forecasts 

of  the  surface  pressure  field  were  computed.  The  forecast 
pressure  change  pattern  was  then  subjectively  compared  with 
the  actual  3  hour  pressure  change  ( isallobaric)  fields  which 
were  analyzed  from  surface  maps.  The  terrain-induced 
vertical  velocity  fields  were  also  subjectively  compared 
with  the  observed  low  level  cloud  patterns.  The  assumption 
is  made  that  the  presence  of  low  cloud  indicates  the 
presence  of  moderately  strong  ascent.  Therefore,  a 
qualitative  comparison  of  the  areal  extent  of  observed  cloud 
patterns  and  vertical  velocity  fields  as  computed  at  the 
lower  boundary  of  the  model  is  an  indication  that  the  lower 
boundary  condition  is  performing  reasonably  well. 
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In  the  preliminary  testing,  it  was  found  that  the 
turning  angles  had  only  a  small  apparent  effect  in  modifying 
the  vertical  velocity  field  and  hence  the  resultant  height 
tendency  fields.  It  was  therefore  decided  to  use  a  fixed 
set  of  turning  angles  in  all  subsequent  testing.  The  choice 
of  these  angles  was  based  partially  on  the  findings  of 
Sawyer  (1959),  and  Greystone  (1962),  and  on  the  summary 
provided  in  Petterssen  (1956),  and  partially  on  the 
subjective  comparisons  of  the  terrain  induced  vertical 
velocity  fields  with  the  observed  low  level  cloud  patterns. 
The  angles  chosen  were  10°  ,20°  and  40°  (to  be  abbreviated 
in  further  discussions  as  o4  :10,20,40)  over  sea,  land  and 
mountains  respectively.  An  exhaustive  testing  of  the  effect 
of  variation  of  the  turning  angles,  employing  many  12  hour 
forecasts,  was  beyond  the  scope  of  this  work. 

The  reduction  factor  r  was  also  tested  in  the  same 
manner  as  c<  and  a  great  variation  in  the  accuracy  of  the 
forecast  was  observed  when  r  was  varied.  As  a  result  of 
this  testing,  the  best  reduction  factors  for  use  in 
determining  the  boundary  layer  wind  were  1.0  over  the  sea, 
0.9  over  land  and  0.8  over  mountains  (abbreviated 
r:1.0,.9,.8).  Since  the  f rictionally-induced  vertical 
velocity  is  proportional  to  the  square  of  the  boundary  layer 
wind,  and  hence  to  the  square  of  the  reduction  factors,  the 
vertical  velocity  is  very  sensitive  to  small  changes  in  r. 
Limitations  on  available  computer  time  did  not  permit  the 
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complete  testing  of  the  reduction  factor  in  1 2  hour 
forecasts.  A  few  runs  with  different  sets  of  r  are  compared 
in  the  following  sections  to  illustrate  the  effect  of  the 
reduction  factor  on  the  predicted  flow  patterns. 

Using  the  above  set  of  turning  angles  and  reduction 
factors,  other  variations  in  the  lower  boundary  condition 
were  then  tested.  Twelve  hour  forecasts  were  run  to  test 
the  effects  of  smoothing  of  the  terrain,  the  use  of 
extrapolated  versus  850  mb  geostrophic  winds  in  calculation 
of  ,  and  the  level  at  which  originates,  either 

terrain  height  or  1000  mb.  The  results  of  these  tests  are 
illustrated  in  Section  5.3. 

In  order  to  obtain  a  quantitative  comparison  of  the 
forecasts,  root  mean  square  (RMS)  errors  between  the 
forecast  and  observed  (verification)  geopotential  height 
fields  were  calculated  at  the  four  primary  forecast  levels 
in  the  model.  The  standard  formula  for  this  quantity 
(Panofsky  and  Brier,  1968)  is 


(5.1) 


where  F  is  the  forecast  geopotential  height  at  point  k  and 
0  is  the  observed  value.  This  statistic  was  calculated  for 
the  entire  grid  and  on  the  subsection  of  the  grid  centered 
over  Alberta  as  indicated  in  Figure  1.  Over  the  whole  grid, 
the  number  of  points  was  483  while  over  the  Alberta  section 
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it  was  35.  The  forecast  geopotential  height  and  surface 
pressure  field  changes  at  the  grid  point  closest  to  Edmonton 
were  also  compared  with  the  12  hour  changes  in  those  fields 
observed  at  Edmonton.  In  the  following  discussions  of  the 
different  runs  made  for  Case  1 ,  these  statistical  measures 
will  be  compared  along  with  the  forecast  height  and  pressure 
patterns. 

5.3  Case  1 :  January  5 »  1972 

Figure  8  shows  the  objective  analysis  of  the  850  and 
500  mb  geopotential  height  fields  at  the  initial  time  in 
Case  1.  An  850  mb  low  and  500  mb  trough  were  initially 
moving  eastwards  across  the  Gulf  of  Alaska  toward  Whitehorse 
(XY)  .  A  trough  at  850  mb  and  also  at  700  mb  (not  shown) 
extended  south  from  the  low  to  just  west  of  Vancouver 
Island.  A  large  area  of  warm  air  advection  extended  ahead 
of  the  low  into  the  southern  Yukon,  central  British  Columbia 
and  into  western  Alberta. 

In  Figure  9,  the  initial  objective  analysis  of  the 
surface  pressure  field  is  compared  with  the  surface  pressure 
field  which  was  extrapolated  from  the  objective  analyses 
aloft.  The  actual  analysis  (top)  indicated  that  a  deep  (968 
mb)  surface  low  was  initially  present  slightly  south  of  the 
position  of  the  low  at  the  850  mb  level.  In  the 
extrapolated  field,  this  feature  is  shifted  northward  but  it 
is  of  approximately  the  same  depth  as  the  analyzed  cyclone. 
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Figure  8.  Initial  isobaric  analysis  at  850  and  500  mb  in 
Case  1.  Geopotential  heights  are  in  decameters.  Contour 
interval  is  6  decameters. 
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Figure  9.  Initial  surface  analysis  and  the  initial  surface 
pressure  field  extrapolated  to  mean  sea  level  in  Case  1. 
Pressures  are  in  mb.  Contour  interval  is  4  mb. 
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The  lee  trough  extending  into  Alberta  is  in  much  the  same 
postion  as  that  of  the  actual  analysis.  The  trough  off  the 
British  Columbia  coast,  however  is  a  much  more  significant 
feature  in  the  extrapolated  analysis,  probably  because  it 
reflects  the  trough  analyzed  at  higher  levels.  The 
overemphasis  of  this  trough  in  the  extrapolated  analysis 
indicates  that  perhaps  the  bogussing  aloft  over  the  Eastern 
Pacific  may  not  have  been  quite  correct.  In  spite  of  the 
obvious  differences  in  the  two  surface  pressure  fields  in 
the  data  sparse  areas ,  the  extrapolation  technique  appears 
to  give  a  reasonable  representation  of  the  mean  sea  level 
pressure  field  at  the  initial  time. 

The  RMS  error  between  the  analyzed  and  extrapolated 
pressure  fields  was  calculated  using  equation  (5.1).  Over 
Alberta,  this  error  was  3.6  mb  while  over  the  whole  map,  the 
error  increased  to  5.0  mb.  These  errors  are  in  part  due  to 
the  bogussing  problems  discussed  in  Chapter  III,  and  in  part 
due  to  errors  inherent  in  the  extrapolation  technique.  An 
estimate  of  the  magnitude  of  the  latter  errors  is  given  in 
Appendix  A. 

The  analyses  of  the  geopotential  and  surface  pressure 
fields  for  0000  GMT  January  6,  1972  (the  verification  fields 
for  the  12  hour  forecasts)  are  shown  in  igures  10  and  11. 
The  850  mb  analysis  indicates  that  a  significant  change  in 
the  low  level  circulation  pattern  had  occured  over  Alberta. 

A  sharp  trough  has  been  formed,  the  axis  of  which  lies 
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Figure  10.  Verification  analysis  at  850  and  500  mb  in  Case 
1.  Heights  in  decameters. 
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through  north  central  Alberta.  In  the  12  hour  period,  the 
850  mb  height  fall  at  Edmonton  (EG)  was  100  gpra.  The  closed 
low  initially  in  the  Gulf  of  Alaska  had  moved  slowly  inland 
and  filled  by  70  gpra  while  the  trough  through  central 
Manitoba  had  shifted  eastward.  At  500  mb,  the  flow  over 
Alberta  shifted  from  northwesterly  to  west-north westerly 
with  the  approach  and  broadening  of  the  ridge  initially 
present  through  central  British  Columbia.  The  trough 
initially  in  the  Gulf  of  Alaska  is  no  longer  present  and  a 
strong,  essentially  zonal  flow  is  evident  all  across  the 
west»central  portion  of  the  grid. 

The  extrapolated  surface  pressure  field  is  presented  in 
Figure  11  along  with  the  surface  analysis  which  is  a 
subjective  hand  analysis  of  the  observed  surface  pressure 
field.  No  objective  analysis  of  this  surface  pressure  field 
was  attempted.  Again  it  is  evident  that  the  extrapolation 
technique  performs  well  with  a  twin-centered  surface  low  in 
Alberta  lying  only  slightly  to  the  northwest  of  its  correct 
position.  The  depth  of  this  feature  has  been  correctly 
extrapolated,  as  has  the  strong  pressure  gradient  lying 
along  the  Bocky  Mountains  from  west-central  Alberta  to 
southern  Montana.  The  extrapolated  analysis  also  indicates 
the  presence  of  a  secondary  low  south  of  Vancouver  Island. 
This  weak  feature  is  apparently  associated  with  the 
occluding  frontal  wave  on  the  actual  surface  analysis,  the 
apex  of  which  was  analyzed  to  be  in  that  area. 
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Figure  11.  Verification  analysis  at  the  surface  and  the 
surface  pressure  field  extrapolated  to  mean  sea  level  from 
the  upper  level  verification  fields.  The  contour  interval 
is  8  mb. 
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5. 4  Comparison  of  Several  Forecast  Runs  Case  \ 

5.4.1  Smoothing  of  the  Topography:  Runs  1A,  IB,  1C 

The  topographical  field  used  in  this  study  was  that 
extracted  by  Schram  (1974)  from  aerological  charts.  While  a 
great  deal  of  smoothing  is  inherent  in  the  extraction 
technique  used,  the  topography  is  quite  steep  in  many  places 
and  slopes  as  large  as  8x10”3  (dimensionless)  exist.  From  a 
scale  analysis  of  quasi-geostrophic  motions,  Haltiner  (1971) 
indicates  that  the  maximum  permissible  slope  should  be  of 
the  order  of  1x10~3.  On  the  other  hand,  a  similar 
calculation  by  Sawyer  (1959)  of  the  ratio  of  vertical  to 
horizontal  dimensions  in  orographically  induced 
disturbances,  suggested  that  slopes  as  large  as  8x10~3  are 
reasonable.  In  the  model  used  by  Greystone  (1962) ,  the 
maximum  slopes  of  his  smoothed  terrain  were  as  large  as 
7x1 0~3.  With  these  conflicting  ideas,  it  then  seems 
desirable  to  test  the  sensitivity  of  the  model  to  the 
steepness  of  the  underlying  terrain. 

In  this  section,  a  comparison  is  made  of  the  influence 
of  two  different  topographical  fields  on  the  forecast.  One 
forecast  uses  the  original  topography,  (Run  1  A)  ,  and  a 
second  (Run  IB)  uses  a  topographical  field  which  is  obtained 
by  means  of  a  single  application  of  the  spatial  smoothing 
operator  described  in  Appendix  B,  to  the  original 
topography.  In  a  third  run,  (1C) ,  the  boundary  layer 
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vertical  velocity  was  calculated  using  the  smoothed 
topography  but  originates  at  1000  mb.  In  Runs  1A  and 
IB,  originates  at  the  particular  topographical  surface 
used.  These  three  runs  demonstrate  the  importance  of  the 
topography  in  the  calculation  of  the  divergence  fields  at 
the  lowest  levels  in  the  model. 

Figure  12  shows  the  smoothed  terrain  field  (converted 
to  pressure  in  the  U.S.  Standard  Atmosphere)  in  which  the 
maximum  terrain  slope  has  been  reduced  to  about  2.5x10-3.  A 
comparison  of  Figure  12  with  the  original  topographical 
field  of  Figure  4  shows  that  the  heights  of  some  ‘peaks'  in 
mountainous  areas  were  reduced  by  up  to  50  mb  by  the 
smoothing  operator. 

In  Runs  1A  to  1C,  the  boundary  layer  winds  were  all 
calculated  through  extrapolation  as  described  in  Section 
4.1.  The  turning  angles  : 10,20,40  and  reduction  factors 
r: 1.0,o 9,. 8  were  used  to  calculate  the  frictional  component 
of  Ccij  .  The  peak  values  of  the  total  vertical  velocity  at 
the  lower  boundary  at  the  initial  time  were  reduced  by  up  to 
20%  when  the  smoothed  boundary  was  used,  but  the  shapes  of 
the  fields  were  not  noticably  altered. 

Quite  significant  differences  in  the  forecast 
geopotential  height  and  surface  pressure  fields  are  apparent 
when  Runs  1A  and  IB  are  compared  in  Figures  1 3  to  15.  A 
comparison  of  the  350  mb  forecast  fields  in  Figure  13 
indicates  that  both  runs  tend  to  fill  the  Gulf  of  Alaska  low 


' 


' 


bT0 


65 


Figure  12,  Smoothed  terrain  pressure  field  labelled  in  mb. 
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Figure  13.  12  hour  forecast  (PROG)  at  850  mb  for  Runs  1A 

and  IB.  Heights  are  in  decameters. 
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and  generate  a  lee  trough  at  850  mb  through  east-central 
Alberta.  However  in  both  cases,  these  tendencies  are 
underforecast  in  12  hours  (C.f.  Figure  10).  A  significant 
difference  in  the  RMS  errors  for  the  850  mb  forecasts  over 
Alberta  is  indicated  in  Table  3  and  the  most  accurate 
forecast  is  obtained  with  the  unsmoothed  terrain.  Figure  14 


CASE  1 

RMS  ERRORS 

MAP/ALBERTA 

FIELD  CHANGES 

FORECAST/ACTUAL 

RUM 

850  mb 

50  0  mb 

850  mb 

5  00  mb 

HSL 

(gpm) 

(gpm) 

(gpm) 

(gpm) 

(mb) 

1  A 

20/1  6 

31/58 

-83/- 100 

80/35  - 

2 6/- 14 

1  B 

19/36 

32/53 

-44/- 100 

52/35 

-9/- 14 

1  C 

25/58 

29/47 

-36/- 100 

56/35 

-6/-1 4 

Table  3.  Root- mean- square  (RMS)  errors  for  Run  1 
calculated  for  the  whole  grid  (MAP)  and  over  the  35  point 
Alberta  subsection  (ALBERTA) .  Field  changes  at  grid  point 
12,14  are  compared  with  the  observed  changes  at  Edmonton, 
Alberta. 

compares  the  500  mb  forecasts  of  Runs  1A  and  IB.  While  the 
differences  between  the  two  are  small,  there  is  a  tendency 
for  both  forecasts  to  overintensify  the  500  mb  trough  as  it 
moves  inland  from  the  Gulf  of  Alaska.  As  a  result  of  this 
spurious  intensification,  the  upper-level  height  gradient 
over  western  Alberta  was  forecast  to  be  somewhat  in  excess 
of  its  actual  value.  The  incorrectly  forecast 
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Figure  14.  12  hour  forecast  at  500  mb  for  Runs  1A 

Heights  are  in  decameters. 


and  IB. 
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intensification  in  the  middle  to  high  troposphere  was 
observed  in  all  runs  performed  with  Case  1 .  The  sources  of 
this  error  will  be  discussed  below. 

At  the  surface,  both  Runs  1A  and  IB  forecast  the 
development  of  a  lee  cyclone  in  north-central  Alberta,  as 
illustrated  in  Figure  15.  This  feature  is  much  better 
defined  in  Run  1A,  and  the  forecast  position  of  the  low 
center  was  within  about  one  grid  length  (200  kra)  of  the 
actual  position  indicated  in  Figure  11.  The  major 
difference  between  the  two  runs  was  in  the  pressure  change 
forecast  in  the  lee  of  the  mountains.  This  difference  is 
illustrated  in  Table  3,  where  the  12  hour  pressure  change 
forecast  with  the  original  terrain  (1A)  was  almost  double 
the  observed  value,  while  with  the  smoothed  topography, 

(IB),  the  pressure  change  was  about  one-half  of  that 
observed. 

From  a  comparison  of  Run  1A  and  Run  IB,  it  is  apparent 
that  the  smoothness  of  the  topographic  surface  has  a 
significant  effect  on  the  formation  and  rate  of 
intensification  of  the  lee  cyclone  in  this  numerical  model. 
The  reason  for  this  effect  becomes  apparent  when  the  nature 
of  the  divergence  term  at  850  mb  in  the  vorticity  equation 
is  examined.  The  finite-difference  analogue  for  this  term 
at  850  mb  is  approximately  "f  ^ ^ —  and  smoothing  of 

the  terrain  will  result  in  a  decrease  of  Cv)<j  and  an 

.  Thus,  this  terra  will  generally  tend  to 


increase  m 
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Figure  15.  12  hour  extrapolated  surface  pressure  forecasts 

for  Runs  1A  and  IB.  Pressures  are  in  mb  and  the  contour 
interval  is  8  mb. 
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decrease  in  the  lee  of  the  mountains  due  to  both  effects  and 
vorticity  production  will  be  reduced.  If  p,  is  set  equal 
to  1000  mb,  i.e.  £0^  originates  at  the  1000  mb  surface, 
the  vorticity  production  is  further  reduced. 

The  twelve  hour  forecast  850  mb  and  surface  pressure 
fields  for  Run  1C  are  shown  in  Figure  16.  The  most  apparent 
features  of  the  850  mb  field  are  the  underprediction  of  the 
deepening  of  the  Alberta  trough  and  the  small  amount  of 
filling  of  the  Gulf  of  Alaska  low  as  it  moved  inland.  The 
surface  pressure  field  shows  no  formation  of  a  closed 
cyclone  in  Alberta  and  the  slightly  intensified  lee  trough 
was  shifted  to  the  north  of  the  position  predicted  in  Runs 
1 A  and  IB.  The  statistics  in  Table  3  verify  the  inaccuracy 
of  this  forecast  at  low  levels  over  Alberta.  It  is  of 
interest,  however,  to  note  that  the  RMS  error  at  500  mb  was 
lower  in  Run  1C  than  in  either  of  the  other  two  runs. 

Cressman  (1963)  observed  only  small  changes  in  the 
hemispheric  prediction  pattern  when  comparisons  were  made  of 
500  mb  forecasts,  with  at  1000  mb  or  at  jPL  .  A 

similar  trend  is  noted  in  the  current  study  at  500  mb  (Table 

3)  but  at  850  mb,  and  at  the  surface,  the  forecast  is  quite 

inferior  when  the  crude  approximation  of  6-tj  at  1000  mb  is 

employed.  In  the  cases  tested  by  Danard  (1966)  using  a 

similar  multilevel  quasi- geostrophic  model,  poor  forecasts 
were  obtained  in  mountainous  areas;  a  result  that  is 
apparently  due  in  part  to  the  fact  that  in  his  model,  63a 
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Figure  16.  12  hour  forecast  at  850  mb  and  the  extrapolated 

surface  forecast  for  Run  1C. 
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originated  at  1000  mb. 

5.4.2  Friction  and  Topography:  Runs  2A,  2B  and  2C 

In  order  to  determine  the  relative  importance  of 
frictional  convergence  and  orographic  lift  in  the  12  hour 
forecasts,  three  runs  were  made.  Run  2A  was  a  forecast  with 
a  completely  smooth  and  flat  lower  boundary.  The  boundary 
layer  vertical  velocity  CO^  was  set  to  zero  everywhere  over 
the  grid,  and  the  lower  boundary  pressure  was  set  to  1000 
mb.  Run  2B  illustrates  the  effect  of  the  orographic 
component  of  vertical  velocity  only.  CO  ^  Mas  calculated 
using  an  extrapolated  horizontal  wind  and  the  smoothed 
topography  described  in  the  previous  section.  In  Run  2C, 
the  full  frictional  effect  was  included  in  addition  to  the 
orographic  influence.  The  frictional  vertical  velocity 
was  calculated  using  the  extrapolated  wind  in  its  unreduced 
form  (r : 1 . 0, 1 . 0, 1 . 0)  and  the  standard  set  of  turning  angles 
oC  :1 0,20, 40  as  well  as  the  smoothed  topography  used  in  Run 
2B. 


In  all  three  of  these  runs,  the  vorticity  field  was 
explicitly  smoothed  every  two  iterations  using  the  smoothing 
operator  described  in  Appendix  B,  as  compared  to  Funs  1A,  IB 
and  1C  where  the  smoothing  was  done  every  12  iterations. 

This  was  done  to  see  whether  more  frequent  smoothing 
significantly  altered  the  forecast  fields  or  RMS  errors. 

The  resulting  forecast  fields  are  shown  in  Figures  17 
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through  20,  and  the  statistical  computations  are  summarized 


RMS  ERRORS  FIELD  CHANGES 

CASE  1  MAP/ALBERTA  FORECAST/ACTUAL 


RON 

850  mb 
(9  Pm) 

500  mb 
(gpm) 

850  mb 
(gpm) 

5  00  mb 
(gpm) 

MS  L 

(mb) 

2  A 

42/51 

45/73 

-63/- 1 00 

53/35 

-9/-1 4 

2  B 

36/59 

39/70 

-28/-100 

51/35 

-5/-14 

2  C 

2  0/4  4 

26/47 

-37/- 100 

60/35 

-7/-1 4 

Table  4.  R  MS  errors  and  field  changes  at  Edmonton  for 
Run  2 . 

in  Table  4.  Discussion  will  again  center  on  the  low  level 
development,  and  a  brief  discussion  of  some  of  the  initial 
CO  fields  will  be  presented.  For  the  sake  of  brevity,  a 
comparison  of  only  the  500  mb  statistics  will  be  attempted, 
the  forecast  flow  patterns  at  this  level  being  visually 
little  different  from  those  of  Runs  1A  and  IB. 

The  850  mb  and  surface  pressure  forecasts  for  the  flat 
frictionless  boundary  condition  of  Run  2A  are  illustrated  in 
Figure  17.  As  is  to  be  expected  with  such  a  poor  boundary 
condition,  no  lee  cyclogenesis  is  predic_ed  and  the  steady 
deepening  of  the  850  mb  and  surface  low  centers  is  apparent. 
The  inland  propagation  of  the  850  mb  and  surface  troughs  is 
clearly  indicated.  Both  move  eastward  a  distance  of  1000  km 
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Figure  17. 
for  a  flat 


12  hour  forecast  at  850  mb  and  at  the  surface 
rictionless  boundary  condition. 
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in  the  12  hour  period  while  the  low  in  the  Gulf  of  Alaska 
moves  very  slowly  eastward.  The  forecast  surface  pressure 
change  at  Edmonton  reflects  the  approach  of  the  trough  at 
the  surface  and  aloft.  The  RMS  errors  for  this  run  are 
large  and  further  indicates  that  in  Case  1 ,  the  lower 
boundary  has  a  major  influence  in  the  development  of  the  low 
level  features. 

Figure  18  shows  the  orographic  component  of  the 
terrain~induced  vertical  velocity  field  at  the  initial  time 
in  Run  2B,  as  well  as  the  CO  field  at  600  mb.  At  the  lower 
boundary,  the  orographic  lift  is  strongest  along  the  west 
coast  of  British  Columbia,  while  several  areas  of  subsidence 
are  evident  to  the  lee  of  the  Rocky  Mountains  through 
northeastern  British  Columbia,  southern  Alberta  and  on  the 
Hontana~*Wyoming  border.  The  maximum  ascent  computed  at  the 
southern  tip  of  the  Alaskan  Panhandle  was  about  12  cm  s~l. 
The  subsidence  in  northeastern  British  Columbia  had  a 
maximum  value  near  3.5  cm  s~4  while  a  more  intense  area  of 
subsidence  associated  with  the  flow  off  the  south  coast  of 
Alaska  had  a  peak  value  of  nearly  8  cm  s”1. 

At  600  mb,  the  CO  pattern  is  typical  of  that  observed 
with  all  combinations  of  the  lower  boundary  condition  at  the 
initial  time  in  the  Case  1  forecasts.  The  maximum  area  of 
ascent  over  the  Alaskan  Panhandle  is  associated  with  the  500 
mb  trough  in  the  Gulf  of  Alaska.  Although  not  indicated  in 
Figure  18,  an  area  of  very  weak  subsidence  was  initially 
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Figure  18.  Orographically  induced  omega  field  at  the 
surface  (bottom)  and  the  omega  field  at  600  mb  (top)  for  the 
initial  time  in  Run  2B.  At  the  surface,  w=~0.9co  where  w 
is  the  vertical  velocity  in  cm  s~ 1  and  oo  is  the  'vertical 
velocity'  in  raicrobars  s~ 1 .  At  600  mb,  w=~1.27co  . 

Negative  values  of  co  denote  ascent  while  positive  values 
denote  subsidence. 
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Figure  19.  12  hour  forecast  at  850  mb  and  at  the  surface 

with  the  orographic  component  of  only  (Run  2B)  . 
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present  over  Alberta  at  this  time. 

In  Figure  19,  the  low  level  forecasts  of  Run  2B  are 
shown,  and  there  is  some  indication  of  an  improvement  over 
the  fields  forecast  in  Run  2A .  The  850  mb  trough  which  was 
observed  to  be  moving  through  British  Columbia  in  Run  2B  is 
now  considerably  weakened  and  a  lee  trough  is  generated  in 
northern  Alberta.  The  rate  of  deepening  of  the  850  mb  low 
is  reduced,  but  the  tendency  is  still  of  the  wrong  sign. 

The  forecast  surface  pattern  shows  some  evidence  of  the 
generation  of  a  weak  cyclone  in  northern  Alberta,  but  the 
surface  cyclone  center  in  the  southwestern  Yukon  is  the 
dominant  feature.  Strong  anticyclogenesis  is  also  occuring 
in  eastern  Idaho.  It  is  associated  with  the  extensive  area 
of  strong  ascent  induced  in  that  region  at  the  lower 
boundary  of  the  model  (see  Figure  18).  This  ascent  results 
in  the  calculation  of  a  large  negative  vorticity  tendency  at 
850  mb  through  the  term  -f  . 

When  the  frictional  component  of  vertical  velocity  is 
included  in  Run  2C,  a  much  more  reasonable  forecast  is 
achieved.  Figure  20  shows  the  full  vertical  velocity 
+  COt  at  the  initial  time.  Several  features  are  of 
particular  importance.  The  term  6J ^  is  quite  large,  as  is 
evident  from  a  comparison  of  Figures  18  s.nd  20.  Along  the 
west  coast  of  British  Columbia,  CO^  is  of  the  same  sign  as 
CO. 'c  and  very  strong  ascent  is  induced  in  this  area. 
Particularly  noteworthy  is  the  strong  band  of  ascent 


. 


*  - 


80 


Figure  20.  Orographically  and  f rictionally-induced  cOe* 
field  at  the  surface  (bottom)  and  the  CO  field  induced  at 
600  mb  (top)  for  the  initial  time  in  Hun  2C. 
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extending  into  the  Pacific  Ocean  south  of  ship  4YP.  This 
band  is  associated  with  the  surface  trough  and  frontal 
system  indicated  in  Figure  9.  The  vertical  velocities  are 
quite  large  within  this  band,  attaining  peak  values  of  about 
11  cm  s~l.  West  of  Edmonton,  subsidence  has  been  enhanced 
by  the  inclusion  of  friction.  In  many  other  areas,  <LO.£  is 
larger  than  and  opposite  in  sign  to  GOt  .  This  is  most 
apparent  in  central  Idaho  where  the  orographic  lift 
indicated  in  Figure  18  has  been  masked  by  the  large 
f rictionally-induced  subsidence. 

At  600  mb,  the  shape  of  the  CO  pattern  at  the  initial 
time  is  quite  similar  to  that  of  Figure  18.  A  notable 
increase  in  the  middle  level  CO  field  is  evident  however. 

At  the  center  of  the  ascent  area  over  Whitehorse  (XY) ,  the 
peak  value  was  increased  by  about  30%  owing  to  the  enhanced 
low  level  ascent. 

The  forecast  850  mb  and  surface  pressure  fields  for  Run 
2C  are  shown  in  Figure  21  where  the  effects  of  the 
frictional  convergence  are  clearly  indicated.  At  850  mb  the 
low  in  the  Gulf  of  Alaska  has  been  completely  filled.  In 
spite  of  the  nearly  correct  height  tendency  forecast  in  this 
area,  the  formation  of  the  lee  trough  in  Alberta  is  still 
not  very  well  handled.  The  forecast  850  mb  height  change  at 
Edmonton  was  much  less  than  the  observed  value  (see  Table 
4)  .  The  RMS  errors  over  Alberta  were  reduced  from  those  of 
Runs  2A  and  2B  while  over  the  whole  map,  this  reduction  was 


' 


82 


Figure  21.  12  hour  forecast  at  850  mb  and  the  surface  v/ith 

full  friction  and  orography  at  the  lower  boundary. 
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quite  large  at  both  850  and  500  mb. 

The  surface  pressure  pattern  also  shows  some 
improvement.  A  closed  lee  cyclone  is  clearly  indicated 
although  the  depth  and  location  of  this  feature  are  not 
quite  correct.  The  effect  of  the  change  in  sign  of  the 
field  over  Idaho  is  also  quite  evident,  and  the  strong 
anticyclogenesis  which  was  predicted  there  in  Run  2B  has 
been  greatly  reduced.  As  well  the  horizontal  pressure 
gradients  over  the  Pacific  near  ship  4 YP  compare  quite 
favorably  with  those  of  the  actual  analysis  for  the  forecast 
verification  time  (Figure  11)  and  they  are  an  improvement 
over  the  strong  gradients  forecast  in  Runs  2A  and  2B. 

In  the  three  runs  described  in  this  section,  the 
importance  of  friction  in  the  development  of  weather  systems 
near  the  mountains  is  apparent.  In  the  model,  strong 
frict ionaliy-induced  ascent  along  the  west  coast  of  British 
Columbia  combined  with  weaker  ascent  aloft  results  in  the 
generation  of  anticyclonic  vorticity  at  850  mb  and  the 
approaching  cyclone  fills.  There  is  however,  a  detrimental 
effect  at  850  mb  in  the  lee  of  the  mountains.  In  Run  2C, 
the  generation  of  anticyclonic  vorticity  due  to  frictional 
ascent  at  the  lower  boundary  prevents  the  proper  deepening 
of  the  850  mb  trough  and  surface  cyclone  in  northern 
Alberta.  Thus  a  slight  adjustment  to  the  frict ionally~ 
induced  CO  field  is  necessary,  and  this  is  accomplished  by 
modifying  the  manner  in  which  the  horizontal  wind  in  the 
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boundary  layer  is  parametarized. 

5.4.3  Small  Frictional  Vertical  Velocity:  Run  3 

In  this  run,  the  importance  of  the  magnitude  of  the 
frictional  convergence  in  the  numerical  model  is  further 
demonstrated.  In  the  calculation  of  c ,  the  turning 
angles  c<:10,20,40  and  the  reduction  factors  r:.5,.2,.2 
were  used  with  an  extrapolated  horizontal  wind  at  the  lower 
boundary.  In  addition,  the  unsmoothed  terrain  was  used  and 
the  vertical  velocity  at  the  lower  boundary  originated  at 
this  terrain  surface. 

The  forecast  850  mb  and  surface  pressure  fields  for 
this  run  are  shown  in  Figure  22  and  the  statistics  are 
listed  in  Table  5.  These  statistics  are  compared  with  those 
of  Run  1 A  since  the  two  runs  are  similar  in  the  sense  that 
they  both  employ  an  unsmoothed  topographical  field  and 
differ  only  in  the  amount  of  frictional  convergence  at  the 
lower  boundary.  Run  3  is  of  interest  because  it  is  the 
worst  low  level  forecast  of  any  of  the  runs  attempted  and  it 
is  one  of  the  worst  at  upper  levels. 

At  850  mb,  the  low  in  the  Gulf  of  Alaska  has  rapidly 
deepened  in  the  12  hour  period.  While  an  850  mb  trough  was 
formed  over  Alberta,  the  rate  of  intensification  of  this 
feature  was  too  rapid  as  is  evident  in  the  field  change 
forecast  at  Edmonton.  The  surface  pressure  forecast  is  a 
reflection  of  the  inaccuracy  at  850  mb.  The  lee  cyclone  was 
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Figure  22.  12  hour  forecast  at  850  mb  and  at  the  surface 

with  a  small  frictionally-induced  OJ  field  at  the  lower 
boundary. 
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generated  in  approximately  the  correct  position  but  it  was 
forecast  approximately  27  mb  too  deep.  A  similar  error  in 
the  forecast  of  the  pressure  change  at  Edmonton  is  indicated 
in  Table  5.  When  these  statistics  are  compared  with  those 
of  Run  1A,  the  differences  are  quite  significant.  Thus  it 


RMS  ERRORS 


FIELD  CHANGES 


CASE  1  MAP/ALBERTA  FORECAST/ACTUAL 

RUN  850  mb  500  mb  850  mb  500  mb  MSL 

(gpm)  (gpm )  (gpm)  (gpra)  (mb) 

3  47/70  37/66  -124/-100  55/35  -35/-14 

1  A  20/16  31/58  ”83/” 100  80/35  -26/-14 


Table  5.  RMS  errors  and  field  changes  at  Edmonton  for 
Run  3  compared  with  those  of  Run  1A. 

is  apparent  that  f rictionally-induced  ascent  plays  an 
important  role  in  the  rate  of  intensification  of  cyclones. 

If  this  ascent  is  not  large  enough  to  oppose  the  orographic 
subsidence  in  the  lee  of  the  mountains,  which  in  this  run  is 
quite  large  with  the  unsmoothed  terrain,  there  is  generally 
strong  vorticity  production  in  the  lee  of  the  mountains 
resulting  in  the  intense  cyclogenesis  observed. 
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5.4.4  The  850  mb  Wind  in  PBL  Computations:  Run  4 

In  all  of  the  forecasts  described  to  this  point,  the 
horizontal  wind  used  in  the  lower  boundary  layer 
calculations  was  obtained  by  extrapolation  from  the  upper 
level  height  fields.  In  this  section,  the  results  of  an 
experiment  are  described  in  which  the  850  mb  geostrophic 
wind  was  used  in  the  PBL.  In  Run  4,  the  turning  angles  and 
reduction  factors  were  identical  to  those  of  Run  IB,  and 
since  the  smoothed  topography  was  also  employed.  Run  4  will 
be  compared  with  Run  IB  to  illustrate  the  differences  which 
occur  in  the  numerical  prediction  when  two  different 
boundary  layer  winds  are  used. 

Figure  23  shows  the  850  mb  and  surface  pressure  fields 
forecast  in  Run  4.  The  low  at  850  mb  was  moved  inland  in  a 
correct  manner  but  the  rate  of  filling  was  quite  slow.  The 
troughing  over  Alberta  was  well  defined  at  this  level,  but 
it  was  slightly  underforecast  both  in  terms  of  the  magnitude 
of  the  height  changes  at  Edmonton  as  well  as  the  sharpness 
of  the  Alberta  trough.  At  the  surface,  the  lee  cyclone  is 
well" organized  and  it  is  forecast  within  3  or  4  mb  of  its 
observed  depth.  However,  the  position  of  the  center  is 
about  five  hundred  kra  northwest  of  the  correct  position. 

When  the  statistics  are  compared  with  those  of  previous 
runs,  it  is  seen  that  overall.  Run  4  produced  the  best 
forecast  with  one  notable  exception  at  850  mb  over  Alberta. 
At  both  850  and  700  mb  (not  shown)  ,  Run  4  exhibited  the 
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Figure  23.  12  hour  forecast  at  850  mb  and  at  the  surface 

with  the  850  mb  geostrophic  wind  used  in  boundary  layer 
calculations. 


RMS  ERRORS 


FIELD  CHANGES 


CASE  1  HAP/ALBERTA  FORECAST/ACTUAL 


RUN 

85  0  mb 

50  0  mb 

850  mb 

500  mb 

MSL 

(gpm) 

(gpm) 

(gpm) 

(gpm) 

(mb) 

4 

15/2  4 

29/54 

-63/” 100 

55/35 

-12/- 1 4 

1  B 

19/3  6 

32/53 

-44/- 100 

52/35 

”9/- 14 

Table  6.  RMS  errors  and  field  changes  at  Edmonton  for 
Run  4  compared  with  those  of  Run  IB. 

lowest  RMS  errors  over  the  whole  map  while  at  5C0  mb  the 
error  was  the  second  lowest  over  both  Alberta  and  the  whole 
map.  The  best  agreement  with  the  observed  pressure  change 
at  Edmonton  was  also  achieved  with  this  run.  These 
statistics  are  somewhat  surprising  because  it  was  not 
expected  that  an  850  mb  wind  would  provide  a  good 
approximation  to  the  PBL  wind,  except  perhaps  in  the  higher 
mountainous  areas.  However  the  difference  in  the  forecasts 
is  likely  due  to  the  fact  that  the  extrapolation  was  to 
terrain  height  since  the  thickness  of  the  PBL  is  assumed  to 
be  zero  in  this  model.  A  model  formulated  with  a  boundary 
layer  of  some  finite  thickness  might  be  expected  to  yield  a 
better  forecast. 

Since  this  run  appears  to  generate  the  best  set  of 
statistics  and  gives  a  low  level  forecast  which  is 
reasonable,  the  lower  boundary  condition  used  in  Run  4  is 
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tested  in  Cases  2  and  3  and  the  results  are  described  below. 
In  all  the  runs  done  on  Case  1,  systematic  errors  were 
apparent  as  the  model  over-intensified  the  500  mb  trough  and 
tended  to  displace  the  low  level  features  somewhat  to  the 
north  of  the  correct  positions  which  are  indicated  in 
Figures  10  and  11.  Similar  errors  were  noted  in  the 
forecasts  done  on  Cases  2  and  3  and  the  sources  of  these 
errors  are  discussed  in  Section  5.7. 

5 . 5  Case  2  :  March  5X  1972 


The  analysis  of  the  initial  and  final  conditions  for 
this  case  are  shown  in  Figures  24  and  25,  and  the  results  of 
two  forecast  runs  are  shown  in  Figures  26  and  27.  The 
statistics  generated  from  these  forecasts  are  tabulated  in 
Table  7.  Two  runs  were  attempted  for  each  of  Cases  2  and  3 
to  show  the  difference  between  two  distinct  sets  of  lower 
boundary  conditions.  Run  1  is  a  forecast  with  the  full 
lower  boundary  condition  as  tested  in  Run  4  of  Case  1  and 
Run  2  has  the  flat  and  frictionless  (simple)  boundary  of  Run 
2A. 


Figure  24  shows  the  850  and  500  mb  analyses  at  the 
initial  time  in  Case  2  (1200  GMT  March  5,  1972).  The 
prominent  feature  at  850  mb  was  the  closed  low  near  the 
Queen  Charlotte  Islands  which  was  moving  rapidly 
northeastward  at  that  time.  This  feature  was  associated 
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Figure  24.  Analysis  at  the  initial  time  at  850  mb  and  at 
500  mb  for  Case  2. 
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with  a  500  nib  trough  and  low  initially  oriented  along 
longitude  140°  W.  With  the  westerly  flow  present  over  the 
Rocky  Mountains,  there  is  some  indication  of  the  presence  of 
a  lee  trough  at  850  rab,  extending  from  Edmonton  to  east  of 
Great  Falls  (GTF)  and  thence  southeastward  through  central 
Wyoming. 

In  the  following  twelve  hour  period,  the  850  mb  low 
continued  moving  east-northeast  to  lie  on  the  Alberta 
British  Columbia  border  as  shown  in  Figure  25.  The  shape  of 
this  low  and  the  associated  trough  became  elongated  with  the 
major  axis  of  the  trough  extending  from  southern  Alaska 
through  west-central  Alberta,  and  thence  essentially  along 
the  lee  of  the  mountains  through  eastern  Wyoming  and 
Colorado.  The  position  of  a  trough  which  had  extended  south 
of  the  initial  850  mb  low,  remained  fixed  during  the  period 
when  the  low  moved  inland.  At  500  mb,  the  southern  section 
of  the  trough  swung  eastward  through  about  8°  of  longitude 
while  the  northern  portion  of  the  trough  remained  fixed. 

Both  the  850  and  500  mb  lows  deepened  slightly  in  the  12 
hour  period. 

The  12  hour  forecasts  of  the  850  mb  fields  for  the  full 
and  simple  lower  boundary  conditions  are  shown  in  Figure  26. 
In  Run  1,  the  overall  appearance  of  the  orecast  pattern  is 
quite  similar  to  the  actual  verification  field  of  Figure  24. 
The  position  of  the  low  center  was  about  200  km  north  of  its 
observed  position,  but  the  orientation  of  the  trough  was 
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Figure  25.  Analysis  at  verification  time  at  850  mb  and  at 
500  mb. 
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Figure  26.  850  mb  forecasts  with  full  lower  boundary 

condition  (Run  1)  and  flat  frictionless  boundary  condition 
(Run  2)  . 
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quite  similar  to  that  of  Figure  25.  This  is  in  contrast  to 
the  pattern  forecast  in  Run  2  with  the  flat  frictionless 
boundary  condition  discussed  in  Run  2A  of  Case  1.  With  no 
orographic  or  frictional  influence  present,  an  essentially 
circular  low  was  forecast  to  move  inland  at  a  much  slower 
speed  than  that  observed.  As  well,  the  trough  south  of  the 
low  was  forecast  to  move  inland  with  the  low  in  Run  2.  In 
Run  1,  frictional  and  orographic  ascent  combined  to  weaken 
and  slow  this  feature  in  its  passage  through  Washington.  A 
weak  ridge  was  maintained  along  a  line  through  the  Alberta- 
British  Columbia  border  and  south  through  western  Montana. 

The  forecast  depth  of  the  850  mb  low  and  the  intensity 
of  the  circulation  were  quite  different  in  Runs  1  and  2. 
While  the  850  mb  low  was  deepened  in  both  runs,  the  rate  was 
not  as  large  in  Run  1  as  in  Run  2.  The  observed  deepening 
of  the  center  was  only  10  gpm  whereas  in  Run  1  a  fall  of  60 
gpra  was  forecast.  Without  the  presence  of  lower  boundary 
effects  to  provide  a  control  on  the  deepening  rate,  a  change 
in  depth  of  -140  gpm  was  forecast  in  Run  2. 

The  500  mb  forecast  of  Run  1  is  shown  in  Figure  27  and 
the  statistics  pertaining  to  this  case  are  given  in  Table  7. 
The  error  in  the  500  mb  geopotential  pattern  reflects  the 
fact  that  the  speed  of  the  trough  and  5C  i  mb  low  were 
forecast  to  be  too  great.  When  Figures  27  and  25  are 
compared,  it  is  seen  that  the  southern  portion  of  the  major 
500  mb  trough  was  moved  a  distance  of  12°  longitude  in  12 
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Figure  27.  500  mb  forecast  with  full  lower  boundary 

condition „ 


RMS  ERRORS  FIELD  CHANGES 


CASE  2  MAP/ALBERTA 

RON  850  mb  500  mb 

(gpm)  (gpm) 

1  24/43  27/53 

2  40/78  39/74 


FORECAST/ACTUAL 

850  mb  500  mb  MSL 
(gpm)  (gpm)  (mb) 

ryy  ea*  <o»  «*  aa  w  »  b  «•  m  s»  29  ca*  v  »  ce»  «?s»  es*  c»  <s» 

”1 21/- 1 17  -59/- 3 1  -16/-15 
-171/- 117  -64/- 31  -27/“ 1 5 


Table  7«  RMS  errors  and  field  changes  at  Edmonton  for 
the  full  and  simple  lower  boundary  conditions  of  Case  2. 
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hours  or  50%  faster  than  the  actual  speed.  As  well,  the  low 
was  forecast  to  deepen  by  50  gpra  instead  of  the  observed  10 
gpm.  A  similar  error  in  the  500  mb  field  over  Edmonton  is 
indicated  in  Table  7  where  the  forecast  change  was 
approximately  twice  the  observed.  I^hen  the  statistics  for 
Runs  1  and  2  are  compared  at  low  levels,  it  is  seen  that  a 
much  better  forecast  was  achieved  with  the  full  lower 
boundary  condition.  In  fact  the  large  850  mb  height  and 
surface  pressure  changes  were  almost  perfectly  forecast  in 
Run  1.  With  this  lower  boundary  condition,  significant 
skill  is  demonstrated. 

5*6  Case  3_:_  April  27 ,  1  974 

The  synoptic  situation  in  this  case  differs 
significantly  from  that  of  the  other  two  cases,  and  in 
contrast  to  Cases  1  and  2,  only  slight  changes  in  the 
weather  patterns  occured  in  the  12  hour  forecast  test 
period.  The  analyses  at  the  initial  time  (0000  GMT  April 

27,  1974),  are  shown  in  Figure  28.  Over  the  central  portion 
of  the  map,  the  circulation  at  500  mb  was  dominated  by  a 
closed  low,  the  center  of  which  was  moving  slowly  eastward 
across  Oregon  and  Idaho  through  a  distance  of  approximately 
300  km  in  12  hours.  As  this  feature  began  to  cross  the 
Rocky  Mountain  chain  at  the  time  of  the  analysis  of  Figure 

28,  the  center  appeared  to  accelerate  and  by  the  end  of  the 
12  hour  period  it  was  located  500  km  to  the  northeast  as 
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Figure  28.  Geopotential  analysis  at  the  initial  time  at  850 
mb  and  500  mb  for  Case  3. 
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Figure  29.  Geopotential  analysis  at  the  verification  time 
at  850  rab  and  500  rab  for  Case  3. 
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shown  in  Figure  29. 

At  low  levels,  southern  Alberta  was  experiencing  heavy 
snow  and  strong  winds  which  may  be  partially  attributable  to 
the  presence  of  the  trough  at  850  rab  in  that  area.  The  low 
center  shown  at  850  mb  in  Figures  28  and  29  had  formed 
several  hours  prior  to  this  analysis  time,  and  it  did  not 
move  appreciably  in  the  12  hour  period  of  this  study.  With 
the  approach  of  the  500  mb  low  however,  the  weak  trough 
analyzed  to  the  west  of  the  850  mb  low  became  quite  sharp  in 
12  hours.  This  intensification,  combined  with  the  strong 
upslope  gradient  through  southern  Alberta  at  all  levels,  was 
largely  responsible  for  the  severe  weather  experienced. 

The  results  of  two  forecast  runs  are  shown  in  Figures 
30  and  31  and  the  statistics  are  summarized  in  Table  8.  The 
same  sets  of  lower  boundary  conditions  were  tested  in  Case  3 
as  in  Case  2.  In  Figure  30,  the  two  850  mb  forecasts  are 
compared.  Run  1  with  the  full  lower  boundary  condition 
shows  a  slight  filling  by  30  gpm  of  the  850  rab  low.  The 
sharpening  trend  of  the  southern  Alberta  trough  is 
indicated,  and  this  is  slightly  more  pronounced  than  that 
forecast  by  Run  2.  In  Run  2,  the  850  mb  low  was  deepened  by 
30  gpm  and  the  center  of  the  low  was  shifted  westward, 
closer  to  the  500  mb  low  center. 

The  results  of  Run  1  are  shown  in  Figure  31.  At  500  mb 
the  motion  of  the  low  center  was  forecast  well.  Although 
the  speed  of  movement  of  the  low  was  slightly  less  than  that 


n 

■ 


101 


Figure  30.  850  mb  forecasts  with  full  lower  boundary 

condition  (Run  1)  and  the  flat  frictionless  boundary 
condition  (Run  2) . 
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Figure  31.  500  mb  forecast  with  the  full  lower  boundary 

condition. 


RMS  ERRORS  FIELD  CHANGES 


CASE  3  MAP/ALBERTA  FOREC AST/ACTOAL 


RUN 

850  mb 
(gpm) 

50  0  mb 
(gpm) 

850  mb 

(gpm) 

5  00  mb 
(gpm) 

MSL 

(mb) 

1 

14/13 

17/12 

20/19 

-16/-10 

7/5 

2 

19/3  5 

15/22 

-39/19 

-39/- 1 0 

-5/5 

Table  8.  RMS  errors  and  field  changes  at  Edmonton  for 
the  full  and  simple  boundary  conditions  of  Case  3. 
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observed,  the  deepening  of  40  gpm  compares  well  with  the 
observed  deepening  of  30  gpm.  In  Run  2  (not  shown)  the  only 
major  difference  was  that  the  low  was  forecast  to  deepen  by 
over  70  gpm  during  the  same  period. 

A  comparison  of  the  statistics  in  Table  8  shows  much 
the  same  trend  as  that  of  the  other  two  cases.  With  the 
full  lower  boundary  condition,  actual  intensification  is 
somewhat  reduced  at  all  levels.  At  850  mb  over  Alberta, 
this  resulted  in  a  reduction  of  RMS  errors  by  over  60%  when 
compared  to  Run  2.  At  500  mb.  Run  1  again  gave  a  better 
forecast  over  Alberta  than  did  Run  2,  although  over  the 
whole  map,  the  flat  frictionless  lower  boundary  condition 
gave  a  marginally  better  forecast.  The  field  changes  at 
Edmonton  were  handled  extremely  well  with  the  full  lower 
boundary  condition  (Run  1)  but  quite  poorly  in  Run  2  where 
the  wrong  trend  was  forecast  at  low  levels. 

In  these  statistics,  one  of  the  most  apparent  features 
is  the  fact  that  the  RMS  errors  are  generally  much  smaller 
than  those  of  Cases  1  or  2.  This  may  be  partly  due  to  the 
better  analysis  of  the  initial  condition  of  Case  3  provided 
by  the  CMC  objective  analysis.  However  the  major  reason  for 
the  small  error  was  that  only  a  small  12  hour  change  in  the 
weather  systems  occured  and  no  rapid  movement  or 
intensification  of  the  weather  systems  was  observed. 

The  results  of  this  case  verify  the  results  noted  above 
in  the  testing  of  Case  1.  It  is  again  evident  that  the 
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lower  boundary  condition  has  a  major  influence  on  numerical 
prediction.  The  particular  boundary  condition  used  in  Run 
1,  with  850  mb  geostrophic  winds  employed  in  the  boundary 
layer  computations,  is  a  parameterization  which  provides  an 
accurate  12  hour  forecast. 

5 • 7  Systematic  Errors  in  the  Model 

One  of  the  most  serious  errors  in  the  model  was  the 
inaccurate  forecast  of  the  upper  level  isobaric  fields. 

This  error  was  largest  at  500  mb  but  in  Case  1  the  500  mb 
error  was  also  of  some  significance.  While  these  fields 
themselves  were  not  of  interest  in  this  study,  errors  at 
these  levels  appeared  to  be  the  cause  of  systematic  errors 
at  lower  levels.  When  the  300  mb  forecast  fields  were 
analyzed,  the  presence  of  a  small  amplitude  wave  pattern  of 
short  wave  length  (2  or  3  grid  lengths)  was  noted.  This 
pattern  was  most  evident  in  Case  1  but  was  also  present  with 
reduced  amplitude  in  Cases  2  and  3.  There  was  a  very  slight 
indication  of  a  similar  occurrence  at  500  mb  in  all  cases. 
The  amplification  of  the  wave  pattern  with  time  indicates 
that  some  form  of  computational  instability  may  have  been 
present  in  the  model. 

When  RMS  errors  at  300  mb  were  comp  .ted,  they  were 
observed  to  be  much  larger  than  the  error  calculated  at 
lower  levels.  In  the  runs  of  Case  1,  the  RMS  error  over  the 
whole  map  ranged  between  70  and  90  gpm,  while  over  Alberta 
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this  error  increased  to  200  gpm  or  higher  in  many  of  the 
runs.  For  Cases  2  and  3  the  RMS  errors  were  much  smaller  at 
300  mb,  dropping  to  around  80  gpm  over  Alberta  in  Case  2  and 
to  30  to  40  gpm  in  Case  3. 

To  estimate  the  effect  of  large  errors  at  300  mb  in  the 
extrapolated  surface  pressure  field,  the  extrapolation 
program  was  tested  on  real  data.  It  was  observed  that  the 
direct  effect  was  essentially  below  detection  and  that 
changes  in  the  300  mb  height  of  200  gpm  resulted  in  a  change 
in  the  extrapolated  pressure  of  only  0.5  mb.  Thus  even 
though  the  300  mb  field  may  be  in  serious  error,  the  surface 
pressure  field  is  not  directly  affected  by  this. 

Errors  in  the  geopotential  fields  at  lower  levels  are 
induced  by  the  errors  aloft  because  all  levels  of  the  model 
are  related  through  the  omega  field.  Because  of  the  intense 
spurious  wave  pattern  induced  at  300  mb,  large  vertical 
velocities  are  generated  at  400  mb  through  the  solution  of 
the  omega  equation.  Since  this  field  is  influential  in 
baroclinic  development  at  500  mb  through  the  divergence  term 
of  the  vorticity  equation,  the  error  at  300  mb  eventually 
propagates  to  500  mb  and  perhaps  even  lower  with  time.  In 
Case  1  for  example,  the  500  mb  trough  over intensifies  as  it 
moves  inland.  This  occurs  because  of  an  incorrect  forecast 
of  the  high  level  CO  fields,  apparently  caused  by  errors 
at  300  mb.  Since  the  major  region  of  500  mb  error  was  over 
the  southern  Yukon,  there  is  a  tendency  for  the 
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extrapolation  of  an  area  of  artificially  low  pressure  in  the 
northern  portions  of  the  grid.  Hence  the  surface  lee 
cyclone  and  trough  tended  to  be  systematically  shifted 
northward  in  Case  1  ,  in  spite  of  a  quite  accurate  forecast 
at  850  mb. 

There  are  two  possible  reasons  for  the  instability  at 
300  mb.  A  careful  examination  of  the  initial  300  mb  field 
generated  by  the  objective  analysis  program  in  Case  1  showed 
that  a  strong  wind  jet  of  78  m  s~l  was  located  directly 
above  the  500  mb  trough  in  the  Gulf  of  Alaska.  From  a 
consideration  of  the  radiosonde  station  reports  of  the 
300  mb  wind  as  well  as  the  CMC  500  mb  analysis  in  that  area, 
a  wind  of  this  magnitude  appeared  to  be  about  two  times  too 
large.  This  difficulty  with  the  analysis  is  believed  to 
have  been  caused  by  the  inability  to  provide  good  300  mb 
bogus  data  in  that  region.  With  the  CFL  condition  (equation 
4.10)  barely  satisfied  at  the  initial  time,  the  jet  may 
strengthen  to  the  point  that  the  stability  criterion  is 
violated.  Thus,  analysis  errors  at  the  initial  time, 
causing  instability  to  be  initiated,  may  be  one  reason  for 
the  large  error  at  300  mb  and  thereby  the  contamination  of 
the  forecast  at  lower  levels. 

A  second  reason  for  the  errors  at  300  mb  is  the 
specification  of  the  lateral  boundary  condition  for  the 
vorticity  equation.  As  Charney  et  al  (1950)  originally 
indicated,  a  sufficient  condition  for  the  solution  of  the 
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vorticity  equation  is  that  the  height  field  be  specified  on 
all  lateral  boundaries  and  the  vorticity  field  on  the  inflow 
boundary.  In  the  model  used  in  this  study,  the  vorticity  is 
specified  on  all  lateral  boundaries,  including  the  outflow 
boundary,  by  means  of  the  observed  vorticity  change  (as 
obtained  from  the  observed  height  change)  over  a  12  hour 
period.  This  change,  is  assumed  to  be  linear  with  time  and 
on  the  outflow  boundary,  an  inconsistency  results  which  can 
be  described  in  terms  of  a  mismatch  between  the  vorticity 
being  advected  towards  the  boundary  and  the  vorticity  being 
carried  through  the  boundary  by  the  boundary  condition.  The 
result  of  this  mismatch  is  that  the  vorticity  not  carried 
through  the  boundary  is  reflected  back  into  the  grid  at  the 
outflow  boundary.  This  phenomenon  is  described  by  Shapiro 
and  O'  Brien  (1970)  who  give  an  example  in  which  such  a 
mismatch  causes  the  formation  of  large  amplitude  noise  in 
the  vorticity  field  at  the  outflow  boundary.  In  a  paper  by 
Eennick  (1973),  several  experiments  were  described  in  which 
such  a  reflection  caused  the  generation  of  two  grid  length 
oscillations  in  the  height  and  vorticity  fields. 

Thus  it  appears  that  the  oscillations  observed  in  the 
300  mb  field  in  the  forecasts  of  this  study  may  have  been 
caused  by  an  incorrect  specification  of  the  outflow  boundary 
condition.  At  300  mb  the  vorticity  advection  is  large 
because  of  the  strong  gradients  present  at  that  level  with 
the  result  that  the  reflection  problems  are  most  pronounced 


at  that  level. 
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There  are  several  methods  whereby  these  problems  are 
eliminated,  the  best  being  a  Lagrangian  advection  scheme 
described  by  Shapiro  and  O'  Brien  (1970)  in  which  the 
vorticity  on  the  outflow  boundary  is  specified  to  be  that 
advected  by  the  model.  The  alternative  in  most  models  is 
the  use  of  spatial  smoothing  to  eliminate  the  reflected 
noise.  Although  smoothing  was  used  in  the  current  study,  it 
did  not  appear  to  be  sufficient  to  remove  the  noise 
completely.  The  smoothing  was  applied  to  the  vorticity 
field  rather  than  the  geopotential  field  and  this  also  may 
have  been  a  mistake.  However  further  testing  of  solutions 
to  this  problem  were  beyond  the  scope  of  this  study. 
Certainly  if  this  type  of  model  were  to  be  used  for  extended 
length  forecasts,  an  improvement  in  the  specification  of  the 
lateral  boundary  condition  would  be  necessary.  It  appears 
that  an  advection  scheme  for  the  outflow  boundary  could 
suffice. 
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CHAPTER  VI 


SUMMARY  AND  CONCLUSIONS 


A  four-level  model  of  the  atmosphere  was  formulated,  in 
which  the  quasi- geostrophic  vorticity  and  omega  equations 
were  solved  to  obtain  short  term  prediction  of  atmospheric 
motions  over  Western  Canada.  These  equations,  with  pressure 
as  the  vertical  coordinate,  were  solved  numerically  on  a 
grid  encompassing  Western  Canada,  the  Yukon,  parts  of  the 
Northwest  Territories,  Alaska,  the  Eastern  Pacific  Ocean  and 
the  northwestern  United  States.  The  grid  point  spacing  was 
200  km.  To  ensure  non-linear  computational  stability,  it 
was  found  necessary  to  employ  the  Adams-Bashf or th  time 
stepping  procedure  with  a  step  length  of  one-half  hour. 

Three  synoptic  situations  were  chosen  in  which  some  degree 
of  cyclogenesis  in  the  lee  of  the  Rocky  Mountains  had 
occurred.  Using  an  objective  analysis  scheme,  the  data  were 
analyzed  at  850,  700,  500  and  300  mb  for  input  in  the  model. 
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The  data  were  obtained  from  the  Northern  Hemispheric  Data 
Tabulation  summaries  for  the  first  two  cases,  bogus  data 
being  introduced  into  the  objective  analysis  program  in  the 
data  sparse  regions  of  the  Eastern  Pacific.  In  the  third 
case,  the  initial  data  consisted  of  the  output  of  the 
objective  analysis  program  used  by  the  Canadian 
Meteorological  Center  in  Dorval,  Quebec.  This  data  was 
interpolated  to  the  grid  used  in  this  study  from  the 
hemispheric  grid  data  employed  at  CMC. 

Twelve  hour  forecasts  were  computed  for  each  of  the 
three  cases  in  order  to  test  different  parame terizations  of 
the  lower  boundary  condition.  Extensive  testing  of  several 
types  of  lower  boundary  condition  was  carried  out  on  Case  1 
while  in  Cases  2  and  3,  only  limited  testing  was  done  to 
verify  the  results  observed  in  the  forecasts  of  Case  1.  The 
vertical  velocity  at  the  lower  boundary,  cO^  ,  consisted  of 
orographic  and  frictional  components,  calculated  using  the 
geostrophic  wind  within  the  model.  In  the  calculation  of 
the  frictional  component  of  CO ^  ,  the  geostrophic  wind  was 
modified  to  simulate  the  wind  near  the  bottom  of  the 
planetary  boundary  layer.  For  purposes  of  comparison, 
forecasts  were  also  computed  with  a  flat,  frictionless  lower 
boundary  where  the  vertical  velocity  was  zero. 

Vertical  velocities  in  the  boundary  layer  were 
calculated  using  two  different  idealized  terrain  surfaces 
with  different  amounts  of  smoothing.  The  vertical  velocity 
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was  assumed  to  originate  at  the  terrain  surface  but  for 
comparison,  a  forecast  was  computed  in  which  the  vertical 
velocity  was  assumed  to  originate  at  the  10C0  mb  surface. 

In  order  to  obtain  a  quantitative  assessment  of  the  many 
forecasts,  a  root-mean- square  error  between  the  forecast  and 
the  verification  fields  was  calculated.  As  well,  the 
forecast  changes  of  the  geopotential  height  and  surface 
pressure  fields  near  Edmonton,  Alberta  were  compared  with 
the  observed  changes  at  Edmonton. 

It  was  found  that  the  pressure  level  at  which  the 
boundary  layer  vertical  velocity  is  assumed  to  originate  is 
important  in  the  numerical  model.  Differences  in  the 
smoothness  of  the  terrain  field  resulted  in  slightly 
different  forecasts.  A  more  accurate  result  was  achieved  if 
the  maximum  slope  of  the  terrain  was  reduced  and  small  scale 
irregularities  in  the  terrain  were  removed  by  spatial 
smoothing.  If  the  boundary  layer  vertical  velocity 
originated  at  the  1000  mb  surface,  the  forecast  over  Alberta 
deteriorated  noticably. 

A  forecast  run  with  only  orographically- induced 
vertical  velocity  at  the  lower  boundary  was  shown  to  exhibit 
comparatively  little  skill  while  if  friction  was  introduced, 
a  much  more  accurate  prognosis  was  achieved.  Frictional 
vertical  velocities  had  to  be  quite  large  to  obtain  an 
accurate  forecast.  Thus  it  appears  that  frictional  effects 
are  at  least  as  important  as  large  scale  orographic  lift  and 
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subsidence,  in  governing  the  behavior  of  weather  systems 
near  the  Canadian  Rocky  Mountains.  This  is  in  broad 
agreement  with  the  results  of  Newton  (1956) ,  who  found  that 
friction  is  the  dominant  factor  in  the  first  few  hours  of 
lee  cyclogenesis.  In  the  current  model,  the  absence  of  any 
frictional  or  orographic  effects  results  in  a  poor  forecast. 

It  appeared  that  the  use  of  an  850  mb  geostrophic  wind, 
rather  than  a  wind  extrapolated  to  the  lower  boundary, 
resulted  in  a  slightly  more  accurate  forecast.  It  is 
probable  however,  that  more  extensive  testing  of  the 
extrapolation  technique,  and  improved  parameterization  of 
the  f rictionally-induced  vertical  velocities  should  result 
in  the  attainment  of  yet  more  accurate  forecasts. 

This  study  has  demonstrated  that  some  aspects  of  the 
lee  cyclogenesis  process  can  be  successfully  duplicated  with 
a  multi-level  quasi- geostrophic  model.  In  forecasting  the 
lee  cyclone,  it  is  of  utmost  importance  to  formulate  a  lower 
boundary  condition  in  which  frictional  effects  are  properly 
parameterized.  Middle  to  upper  tropospheric  flow  as  well  as 
that  at  lowest  levels  is  modified  if  boundary  layer  effects 
are  properly  taken  into  account. 
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APPENDIX  A 


THE  AITKEN  POLYNOMIAL 


Given  the  set  of  n+1  points  f(xe)^  ,  £x,  ;  *f  (X,)  j  * 
£'X>7,  -pCx„)j,  the  n'th  order  Aitken  polynomial  which 
interpolates  this  set  is  defined  by: 


m  t> 


x„-x 


DET 


>n 


x^-x 


(A.1) 


Using  this  result,  a  second-order  parabolic  interpolation  of 
the  function  f (x)  is: 


(x*-x)£  -  ( 


(A.  2) 


where 


IF^,  =  x^Co  [(x.-x)f(x.)  -  (Xo-x)-f  Cx,)J  <A-3> 


and 
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ffL  =  [(xi-x)fcx„)-CXo-x)-f(xjJ 


(A. 4) 


The  first  derivative  of 


Kn, 


with  respect  to  X  is  simply: 


!Po,l  =KrX,  [(<»-5f)K./-Cx-x)|^,'-|g1  +  Kj  (A.5) 


where 

and 


K,'=  x^.[-fcx,)--fcx.)] 

fu>2  -  [  f  OO  ”  ^ 


The  second  derivative 


R” 


reduces  to: 


Ip"  _  p  (  ^ 

K,a.  ~  ^  V  X,  -  X, 


(A. 6) 


In  the  numerical  model,  the  CO  field  is  assumed  to  have 
a  parabolic  profile  with  respect  to  pressure  and  the  first 
and  second  derivatives  of  CO  are  obtained  from  equations 
(A.  3)  to  (A.  6)  where  and  CJi  -  f(xt).  Then  given  the 

omega  field  (  C09  ,  CO,  ,  CO A  )  at  three  levels  {  pa  ,  pt  ,  )  ,  the 

first  derivative  of  CO  at  level  p  can  be  shown  to  be  of  the 
finite- difference  form: 

<§Tf3  ~~  Pa.- P,  L  P,  -  J°o  ^6c,»  +  "2- P)  ~  +  P)J 

(A. 7) 

while  the  second  derivative  is  simply: 
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o)  __  {L-  (  CJz.  60  *  60,  -  COo 

/Vfi  V^~P0  >,  -  n  j 


(A. 8) 


It  may  be  readily  verified  that  if  these  derivatives  are 
evaluated  at  pressure  P,  and  if  the  differences  -  p  and 
P,  -  Pa  are  equal,  the  standard  centered-dif f erence  forms  for 
the  first  and  second  derivatives  are  obtained. 


As  well  as  using  the  Aitken  polynomial  to  interpolate 
the  CO  field,  a  third-order  form  \ Poia3  is  used  to  approximate 
the  logarithmic  variation  of  geopotential  height  with 
pressure  in  the  vertical.  This  third-order  polynomial  is 
then  used  to  extrapolate  the  mean  sea  level  pressure  to  zero 
geopotential  height.  For  this  case  the  dependent  variable 
is  the  pressure  and  the  independent  variable  is  the  height 
of  the  four  isobaric  surfaces:  850,  700,  500  and  300  mb 
above  each  point  in  the  horizontal  grid. 

In  order  to  obtain  an  estimate  of  the  error  associated 
with  the  use  of  the  interpolating  polynomial  to  extrapolate 
surface  pressures,  the  technique  was  tested  on  a  randomly 
chosen  set  of  18  radiosonde  ascent  reports.  The 
geopotential  heights  of  the  four  isobaric  surfaces  were  used 
in  the  program  and  the  pressure  at  the  L  ight  above  mean  sea 
level  of  the  radiosonde  station  was  extrapolated.  This 
pressure  was  then  compared  with  the  actual  station  pressure 
reported  at  the  radiosonde  station.  The  deviations  from  the 
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observed  pressure  were  in  the  range  of  zero  to  six  mb  and 
the  standard  deviation  for  the  18  soundings  (at  10  different 
stations)  was  calculated  to  be  2.6  mb.  This  error  is 
assumed  to  be  the  approximate  error  existing  in  the 
extrapolation  of  the  mean  sea  level  pressure  maps  in  the 
discussions  in  Case  1. 
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APPENDIX  B 


SPATIAL  SMOOTHING 


Small  scale  irregularities  are  introduced  into 
meteorological  fields  by  the  numerical  techniques  used  in 
the  objective  analysis  programs.  These  irregularities  are 
often  described  as  noise,  which  is  superimposed  upon  the 
meteorologically  significant  wave  patterns.  Noise  of  wave 
length  two-grid-lengths  or  less  is  a  particular  problem 
since  such  noise  is  aliased  into  the  longer  wave  patterns. 
When  numerical  prediction  of  the  meteorological  waves  is 
attempted,  these  aliased  short  waves  will  often  tend  to 
amplify  with  time,  thereby  distorting  or  obscuring  the 
meteorological  phenomena. 

Even  if  this  noise  is  not  initially  present,  it  may  by 
introduced  into  the  prediction  because  of  the  truncation 
errors  associated  with  the  finite-difference  representation 
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of  the  partial  derivatives  in  the  meteorological  equations. 
Thus  it  is  desirable  to  remove  the  noise  initially  and 
perhaps  periodically  from  the  fields.  This  is  accomplished 
by  means  of  spatial  smoothing. 

The  theory  of  smoothing  techniques  will  not  be 
described  here.  The  reader  is  referred  to  discussions  in 
Haltiner  (1971)  and  Asselin  (1966) .  A  two-dimensional 
smoothing  operator  may  be  written  in  finite-difference  form 
as: 


-  %  *  %  rf„ 

where  -f  is  the  smoothed  value  of  the  field  f  at  point 

is  the  five-point  Laplacian  operator  defined  in 
equation  (4.4)  and  S  is  a  smoothing  parameter  which  may  be 
negative.  In  the  current  study,  S=1/2  was  used  in  all 
smoothing.  It  may  be  shown  that  this  value  of  the  smoothing 
parameter  completely  filters  out  the  two-grid-length  noise 
and  damps  out  waves  shorter  than  two-grid-lengths.  Waves 
longer  then  two-grid-lengths  are  also  damped  somewhat,  but 
this  damping  is  not  very  severe  for  the  meteorologically 
significant  waves  of  the  quasi-geostrophic  system. 

Smoothing  was  applied  to  the  vortic  ty  field  and  was 
done  both  at  the  initial  time  and  at  regularly  spaced 
intervals  in  the  forecast  period.  The  frequency  with  which 
smoothing  was  performed  seemed  to  have  little  effect  on  the 
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magnitude  of  the  RMS  errors  in  the  forecasts. 

The  same  smoothing  operator  was  also  applied  to  the 
pressure  and  density  fields  at  terrain  height  (obtained 
using  the  0,  S.  Standard  Atmosphere) . 


